A hybrid particle-continuum method for hydrodynamics of complex fluids
AA hybrid particle-continuum method for hydrodynamics of complex fluids
Aleksandar Donev, ∗ John B. Bell, Alejandro L. Garcia, and Berni J. Alder Center for Computational Science and Engineering,Lawrence Berkeley National Laboratory, Berkeley, CA, 94720 Department of Physics and Astronomy, San Jose State University, San Jose, California, 95192 Lawrence Livermore National Laboratory,P.O.Box 808, Livermore, CA 94551-9900
A previously-developed hybrid particle-continuum method [
J. B. Bell, A. Garcia and S.A. Williams, SIAM Multiscale Modeling and Simulation, 6:1256-1280, 2008 ] is generalizedto dense fluids and two and three dimensional flows. The scheme couples an explicit fluc-tuating compressible Navier-Stokes solver with the Isotropic Direct Simulation Monte Carlo(DSMC) particle method [
A. Donev and A. L. Garcia and B. J. Alder, ArXiv preprint0908.0510 ]. To achieve bidirectional dynamic coupling between the particle (microscale)and continuum (macroscale) regions, the continuum solver provides state-based boundaryconditions to the particle subdomain, while the particle solver provides flux-based bound-ary conditions for the continuum subdomain. This type of coupling ensures both state andflux continuity across the particle-continuum interface analogous to coupling approaches fordeterministic parabolic partial differential equations; here, when fluctuations are included,a small ( < ∗ Electronic address: [email protected] a r X i v : . [ c ond - m a t . s o f t ] O c t I. INTRODUCTION
With the increased interest in nano- and micro-fluidics, it has become necessary to develop toolsfor hydrodynamic calculations at the atomistic scale [1, 2, 3]. While the Navier-Stokes-Fouriercontinuum equations have been surprisingly successful in modeling microscopic flows [4], there areseveral issues present in microscopic flows that are difficult to account for in models relying on apurely PDE approximation. For example, it is well known that the Navier-Stokes equations failto describe flows in the kinetic regions (large Knudsen number flows) that appear in small-scalegas flows [5]. It is also not a priori obvious how to account for the bidirectional coupling betweenthe flow and embedded micro-geometry or complex boundaries. Furthermore, it is not trivial toinclude thermal fluctuations in Navier-Stokes solvers [6, 7, 8, 9], and in fact, most of the timethe fluctuations are omitted even though they can be important at instabilities [10] or in drivingpolymer dynamics [11, 12]. An alternative is to use particle-based methods, which are explicitand unconditionally stable, robust, and simple to implement. The fluid particles can be directlycoupled to the microgeometry, for example, they can directly interact with the beads of a polymerchain. Fluctuations are naturally present and can be tuned to have the correct spatio-temporalcorrelations.Several particle methods have been described in the literature, such as molecular dynamics(MD) [13], Direct Simulation Monte Carlo (DSMC) [14], dissipative particle dynamics (DPD) [15],and multi-particle collision dynamics (MPCD) [16, 17]. Here we use the Isotropic DSMC (I-DSMC)stochastic particle method described in Ref. [18]. In the I-DSMC method, deterministic interactionsbetween the particles are replaced with stochastic momentum exchange (collisions) between nearbyparticles. The I-DSMC method preserves the essential hydrodynamic properties of expensive MD:local momentum conservation, linear momentum exchange on length scales comparable to theparticle size, and a similar fluctuation spectrum. At the same time, the I-DSMC fluid is ideal andstructureless, and as such is much simpler to couple to a continuum solver.However, even particle methods with coarse-grained dynamics, such as I-DSMC, lack the effi-ciency necessary to study realistic problems because of the very large numbers of particles needed tofill the required computational domain. Most of the computational effort in such a particle methodwould, however, be expended on particles far away from the region of interest, where a descriptionbased on the Navier-Stokes equations is likely to be adequate. Hybrid methods are a natural can-didate to combine the best features of the particle and continuum descriptions. A particle methodcan be used in regions where the continuum description fails or is difficult to implement, such asnear suspended structures or complex boundaries, and a more efficient continuum description canbe used to around the particle domain, as illustrated in Fig. 1. This type of hybrid algorithm fitsin the Multi-Algorithm Refinement (MAR) simulation approach to the modeling and simulationof multiscale problems [19, 20, 21]. MAR combines two or more algorithms, each of which is ap-propriate for a different scale regime. The general idea is to perform detailed calculations usingan accurate but expensive (e.g., particle) algorithm in a small region, and couple this computa-tion to a simpler, less expensive (e.g., continuum) method applied to the rest. The challenge isto ensure that the numerical coupling between the different levels of the algorithm hierarchy, andespecially the coupling of the particle and continuum computations, is self-consistent, stable, andmost importantly, does not adversely impact the underlying physics.
Figure 1: A two-dimensional illustration of the use of the MAR hybrid to study a polymer chain (larger redcircles) suspended in an I-DSMC particle fluid. The region around the chain is filled with particles (smallergreen circles), while the remainder is handled using a fluctuating hydrodynamic solver. The continuum(macro) solver grid is shown (thicker blue lines), along with the (micro) grid used by the particle method(thinner blue lines). The fluctuating velocities in the continuum region are shown as vectors originating fromthe center of the corresponding cell (purple). The interface between the particle and continuum regions ishighlighted (thicker red line).
A crucial feature of our hybrid algorithm is that the continuum solver includes thermal fluctu-ations in the hydrodynamic equations consistent with the particle dynamics, as previously investi-gated in one dimension in Ref. [21]. Thermal fluctuations play an important role in describing thestate of the fluid at microscopic and mesoscopic scales, especially when investigating systems wherethe microscopic fluctuations drive a macroscopic phenomenon such as the evolution of instabilities,or where the thermal fluctuations drive the motion of suspended microscopic objects in complexfluids. Some examples in which spontaneous fluctuations can significantly affect the dynamics in-clude the breakup of droplets in jets [22, 23], Brownian molecular motors [24], Rayleigh-Bernardconvection [25], Kolmogorov flows [26, 27], Rayleigh-Taylor mixing [10], and combustion and ex-plosive detonation [28]. In our algorithm, the continuum solver is a recently-developed three-stageRunge-Kutta integration scheme for the Landau-Lifshitz Navier-Stokes (LLNS) equations of fluc-tuating hydrodynamics in three dimensions [9], although other finite-volume explicit schemes cantrivially be substituted.As summarized in Section II, the proposed hybrid algorithm is based on a fully dynamic bidirec-tional state-flux coupling between the particle and continuum regions. In this coupling scheme thecontinuum method provides state-based boundary conditions to the particle subdomain throughreservoir particles inserted at the boundary of the particle region at every particle time step. Duringeach continuum time step a certain number of particle time steps are taken and the total particleflux through the particle-continuum interface is recorded. This flux is then imposed as a flux-basedboundary condition on the continuum solver, ensuring strict conservation [20, 21]. Section IIIdescribes the technical details of the hybrid algorithm, focusing on components that are distinctfrom those described in Refs. [20, 21]. Notably, the use of the Isotropic DSMC particle methodinstead of the traditional DSMC method requires accounting for the interactions among particlesthat are in different continuum cells.In Section IV we thoroughly test the hybrid scheme in both equilibrium and non-equilibriumsituations, and in both two and three dimensions. In Section IV A we study the continuity ofdensity and temperature across the particle-continuum interface and identify a small mismatch oforder 1 /N , where N is the number of particles per continuum cell, that can be attributed to theuse of fluctuating values instead of means. In Section IV B we compute dynamic structure factorsin periodic (“bulk”) and finite quasi two- and one-dimensional systems and find that the hybridmethod seamlessly propagates thermal fluctuations across the particle-continuum interface.In Section IV C we study the diffusive motion of a large spherical buoyant bead suspended ina bead of I-DSMC particles in three dimensions by placing a mobile particle region around thesuspended bead. This example is of fundamental importance in complex fluids and micro-fluidics,where the motion of suspended objects such as colloidal particles or polymer chains has to besimulated. Fluctuations play a critical role since they are responsible for the diffusive motion of thebead. The velocity-autocorrelation function (VACF) of a diffusing bead has a well-known powerlaw tail of hydrodynamic origin and its integral determines the diffusion coefficient. Therefore,computing the VACF is an excellent test for the ability of the hybrid method to capture theinfluence of hydrodynamics on the macroscopic properties of complex fluids.Finally, in Section IV D we study the slow relaxation toward thermal equilibrium of an adiabaticpiston with the particle region localized around the piston. In the formulation that we consider, thesystem is bounded by adiabatic walls on each end and is divided into two chambers by a mobile andthermally insulating piston that can move without friction. We focus on the case when the initialstate of the system is in mechanical equilibrium but not thermodynamic equilibrium: the pressureson the two sides are equal but the temperatures are not. Here we study the slow equilibration ofthe piston towards the state of thermodynamic equilibrium, which happens because asymmetricfluctuations on the two sides of the piston slowly transfer energy from the hotter to the colderchamber.We access the performance of the hybrid by comparing to purely particle simulations, whichare assumed to be “correct”. Unlike in particle methods, in continuum methods we can triviallyturn off fluctuations by not including stochastic fluxes in the Navier-Stokes equations. By turningoff fluctuations in the continuum region we obtain a deterministic hybrid scheme, to be contrastedwith the stochastic hybrid scheme in which fluctuations in the continuum region are consistent withthose in the particle region. By comparing results between the deterministic and stochastic hybridwe are able to assess the importance of fluctuations. We find that the deterministic hybrid gives thewrong long-time behavior for both the diffusing spherical bead and the adiabatic piston, while thestochastic hybrid correctly reproduces the purely particle runs at a fraction of the computationalcost. These examples demonstrate the need to include thermal fluctuations in the continuumsolvers in hybrid particle-continuum methods. II. BRIEF OVERVIEW OF THE HYBRID METHOD
In this section we briefly introduce the basic concepts behind the hybrid method, delegatingfurther technical details to later sections. Our scheme is based on an Adaptive Mesh and Algo-rithm Refinement (AMAR) methodology developed over the last decade in a series of works inwhich a DSMC particle fluid was first coupled to a deterministic compressible Navier-Stokes solverin three dimensions [20, 29], and then to a stochastic (fluctuating) continuum solver in one dimen-sion [21]. This section presents the additional modifications to the previous algorithms necessaryto replace the traditional DSMC particle method with the isotropic DSMC method [18], and afull three-dimensional dynamic coupling of a complex particle fluid [30] with a robust fluctuatinghydrodynamic solver [9]. These novel techniques are discussed further in Section III.Next, we briefly describe the two components of the hybrid, namely, the particle microscopic model and the continuum macroscopic solver, and then outline the domain decomposition used tocouple the two, including a comparison with other proposed schemes. Both the particle algorithmand the macroscopic solver have already been described in detail in the literature, and furthermore,both can easily be replaced by other methods. Specifically, any variant of DSMC and MPCD canbe used as a particle algorithm, and any explicit finite-volume method can be used as a continuumsolver. For this reason, in this paper we focus on the coupling algorithm.
A. Particle Model
The particle method that we employ is the Isotropic DSMC (I-DSMC) method, a dense fluid generalization of the Direct Simulation Monte Carlo (DSMC) algorithm for rarefied gas flows [14] .The I-DSMC method is described in detail in Ref. [18] and here we only briefly summarize it. It isimportant to note that, like the traditional DSMC fluid, the I-DSMC fluid is an ideal fluid, that is,it has the equation of state (EOS) and structure of an ideal gas; it can be thought of as a viscousideal gas. As we will see shortly, the lack of structure in the I-DSMC fluid significantly simplifiescoupling to a continuum solver while retaining many of the salient features of a dense fluid.In the I-DSMC method, the effect of interatomic forces is replaced by stochastic collisionsbetween nearby particles. The interaction range is controlled via the collision diameter D or,equivalently, the density (hard-sphere volume fraction) φ = πN D / (6 V ), where N is the totalnumber of particles in the simulation volume V . The strength of the interaction is controlledthrough the dimensionless cross-section prefactor χ , which is on the order of unity [18]. Stochasticcollisions are processed at the end of every particle time step of duration ∆ t P , and in-betweencollisions each particle i is streamed advectively with constant velocity υ i = ˙ r i . The spatialdomain of the simulation, typically a rectangular region, is divided into micro cells of length L c (cid:39) D , which are used to efficiently identify all particles that are within distance L c of a givenparticle by searching among the particles in neighboring cells (each cell has 3 d neighboring cells,counting itself, where d is the spatial dimension). At each time step, conservative collisions occurbetween randomly-chosen pairs of particles that are closer than a distance D apart; specifically, Note that by a dense fluid we mean a fluid where the mean free path is small compared to the typical fluidinter-atomic distance. a random amount of momentum and energy is exchanged between the two particles with theprobabilities of the collision outcomes obeying detailed balance. The I-DSMC algorithm can beviewed as a stochastic alternative to deterministic hard-sphere molecular dynamics (MD) [31],where hard spheres of diameter D collide when they touch.In addition to the fluid (I-DSMC) particles there may be a number of additional sphericalparticles, which we refer to as beads , suspended in the I-DSMC fluid. The beads interact witheach other and the fluid particles either deterministically, as hard spheres impermeable to andcolliding with touching particles, or stochastically, as permeable spheres colliding stochastically withoverlapping particles. Note that the sphere radii used for determining the fluid-fluid, fluid-bead,and bead-bead interaction distances may, in general, be different. Combining deterministic hard-sphere collisions with stochastic collisions requires a mixed time-driven and event-driven treatment,as in the Stochastic Event-Driven MD (SEDMD) algorithm developed in Ref. [30]. We will usethe SEDMD algorithm in the examples presented in this paper.We have also developed an alternative purely time-driven fluid-bead coupling in which the fluidis allowed to permeate the beads and all of the particle interactions are stochastic. This leads toa much simpler algorithm that can also easily be parallelized. This distinction between hard andpermeable spheres has analogues in other methods for complex fluids. For example, in Lattice-Boltzmann simulations [32], beads can either be modeled as hard spheres (using a bounce-backcollision rule at the lattice sites on the surface of the bead), or, more efficiently and commonly, aspermeable spheres that let the fluid pass through them but experience a frictional force due to thefluid motion (exerting the opposite force back on the lattice sites they overlap with). B. Continuum Model
At length scales and time scales larger than the molecular ones, the dynamics of the particlefluid can be coarse grained [33, 34] to obtain evolution equations for the slow macroscopic variables.Specifically, we consider the continuum conserved fields U ( r , t ) = ρ j e ∼ = (cid:101) U ( r , t ) = (cid:88) i m i p i e i δ [ r − r i ( t )] = (cid:88) i υ i υ i / m i δ [ r − r i ( t )] , (1)where the conserved variables , namely the densities of mass ρ , momentum j = ρ v , and energy e = (cid:15) ( ρ, T ) + ρv , can be expressed in terms of the primitive variables , mass density ρ , velocity v and temperature T ; here (cid:15) is the internal energy density. Here the symbol ∼ = means that we considera stochastic field U ( r , t ) that approximates the behavior of the true atomistic configuration (cid:101) U ( r , t )over long length and time scales (compared to atomistic scales) in a certain integral average sense;notably, for sufficiently large cells the integral of U ( r , t ) over the cell corresponds to the totalparticle mass, momentum and kinetic energy contained inside the cell.The evolution of the field U ( r , t ) is modeled with the Landau-Lifshitz Navier-Stokes (LLNS)system of stochastic partial differential equations (SPDEs) in d dimensions, given in conservativeform by ∂ t U = − ∇ · [ F ( U ) − Z ( U , r , t )] , (2)where the deterministic flux is taken from the traditional compressible Navier-Stokes-Fourier equa-tions, F ( U ) = ρ v ρ vv T + P I − σ ( e + P ) v − ( σ · v + ξ ) , where P = P ( ρ, T ) is the pressure, the viscous stress tensor is σ = 2 η (cid:104) ( ∇ v + ∇ v T ) − ( ∇ · v ) d I (cid:105) (we have assumed zero bulk viscosity), and the heat flux is ξ = µ ∇ T . As postulated by Landau-Lifshitz [34, 35], the stochastic flux Z = · v + Ξ is composed of the stochastic stress tensor Σ and stochastic heat flux vector Ξ , assumed to bemutually uncorrelated random Gaussian fields with a covariance (cid:10) Σ ( r , t ) Σ (cid:63) ( r (cid:48) , t (cid:48) ) (cid:11) = C Σ δ ( t − t (cid:48) ) δ ( r − r (cid:48) ), where C ( Σ ) ij,kl = 2¯ ηk B T (cid:18) δ ik δ jl + δ il δ jk − d δ ij δ kl (cid:19)(cid:10) Ξ ( r , t ) Ξ (cid:63) ( r (cid:48) , t (cid:48) ) (cid:11) = C Ξ δ ( t − t (cid:48) ) δ ( r − r (cid:48) ), where C ( Ξ ) i,j = 2¯ µk B T δ ij , (3)where overbars denote mean values.As discussed in Ref. [9], the LLNS equations do not quite make sense written as a system ofnonlinear SPDEs, however, they can be linearized to obtain a well-defined linear system whoseequilibrium solutions are Gaussian fields with known covariances. We use a finite-volume dis-cretization, in which space is discretized into N c identical macro cells V j of volume V c , and thevalue U j stored in cell 1 ≤ j ≤ N c is the average of the corresponding variable over the cell U j ( t ) = 1 V c (cid:90) V j U ( r , t ) d r = 1 V c (cid:90) V j (cid:101) U ( r , t ) d r , (4)where (cid:101) U is defined in Eq. (1). Time is discretized with a time step ∆ t C , approximating U ( r , t )pointwise in time with U n = (cid:8) U n , ..., U nN c (cid:9) , U nj ≈ U j ( n ∆ t C ) , where n ≥ n ex of micro time steps,∆ t C = n ex ∆ t P .In addition to the cell averages U nj , the continuum solver needs to store the continuum normalflux F nj,j (cid:48) through each interface I = V j ∩ V j (cid:48) between touching macro cells j and j (cid:48) during a giventime step, U n +1 j = U nj − ∆ tV c (cid:88) j (cid:48) S j,j (cid:48) F nj,j (cid:48) , (5)where S j,j (cid:48) is the surface area of the interface, and F nj,j (cid:48) = − F nj (cid:48) ,j . Here we will absorb thevarious prefactors into a total transport (surface and time integrated flux) through a given macrocell interface I , Φ nI = V − c ∆ tS I F nI , which simply measures the total mass, momentum and energytransported through the surface I during the time interval from time t to time t +∆ t . We arbitrarilyassign one of the two possible orientations (direction of the normal vector) for each cell-cell interface.How the (integrated) fluxes Φ nI are calculated from U nj does not formally matter; all that the hybridmethod uses to advance the solution for one macro time step are U nj , U n +1 j and Φ nI . Therefore, anyexplicit conservative finite-volume method can be substituted trivially. Given this generality, wedo not describe in any detail the numerical method used to integrate the LLNS equations; readerscan consult Ref. [9] for further information. C. Coupling between particle and continuum subdomains
The hybrid method we use is based on domain decomposition,and is inspired by Adaptive MeshRefinement (AMR) methodology for conservation laws [36, 37]. Our coupling scheme closely followspreviously-developed methodology for coupling a traditional DSMC gas to a continuum fluid, firstproposed in the deterministic setting in Ref. [20] and then extended to a fluctuating continuummethod in Ref. [21]. The key new ingredient is the special handling of the collisional momentumand energy transport across cell interfaces, not found in traditional DSMC. For completeness, wedescribe the coupling algorithm in detail, including components already described in the literature.0We split the whole computational domain into particle and continuum subdomains, which com-municate with each other through information near the particle-continuum interface I , assumedhere to be oriented such that the flux Φ I measures the transport of conserved quantities fromthe particle to the continuum regions. In AMR implementations subdomains are usually logicallyrectangular patches ; in our implementation, we simply label each macro cell as either a particlecell or a continuum cell based on whatever criterion is appropriate, without any further restrictionson the shape or number of the resulting subdomains. For complex fluids applications, macro cellsnear beads (suspended solute) and sometimes near complex boundaries will be labeled as particlecells. The continuum solver is completely oblivious to what happens inside the particle subdomainand thus it need not know how to deal with complex moving boundaries and suspended objects.Instead, the continuum solution feels the influence of boundaries and beads through its couplingwith the particle subdomains.The dynamic coupling between particle and continuum subdomains is best viewed as a mu-tual exchange of boundary conditions between the two regions. Broadly speaking, domain-decomposition coupling schemes can be categorized based on the type of boundary conditionseach subdomain specifies for the other [38]. Our scheme is closest to a state-flux coupling schemebased on the classification proposed in Ref. [38] for incompressible solvers (the term “velocity-flux” is used there since velocity is the only state variable). A state-flux coupling scheme is onein which the continuum solver provides to the particle solver the conserved variables U in thecontinuum reservoir macro cells near the particle-continuum interface I , that is, the continuumstate is imposed as a boundary condition on the particle region. The particle solver provides tothe continuum solver the flux Φ I through the interface I , that is, the particle flux is imposed asa boundary condition on the continuum subdomain. This aims to achieve continuity of both statevariables and fluxes across the interface, and ensures strict conservation, thus making the couplingrather robust. Note that state/flux information is only exchanged between the continuum andparticle subdomains every n ex particle (micro) time steps, at the beginning/end of a macro timestep. A more detailed description of the algorithm is given in Section III.
1. Comparison with other coupling schemes
There are several hybrid methods in the literature coupling a particle method, and in partic-ular, molecular dynamics (MD), with a continuum fluid solver [39]. There are two main typesof applications of such hybrids [40]. The first type are problems where the particle description1is localized to a region of space where the continuum description fails, such as, for example, acomplex boundary, flow near a corner, a contact line, a drop pinchoff region, etc.. The second typeare problems where the continuum method needs some transport coefficients, e.g., stress-strainrelations, that are not known a priori and are obtained via localized MD computations. In themajority of existing methods a stationary solution is sought [41], and a deterministic incompress-ible or isothermal formulation of the Navier-Stokes equations is used in the continuum [40, 42]. Bycontrast, we are interested in a fully dynamic bidirectional coupling capable of capturing the fullrange of hydrodynamic effects including sound and energy transport. We also wish to minimizethe size of the particle regions and only localize the particle computations near suspended objects,making it important to minimize the artifacts at the interface. As we will demonstrate in thiswork, including thermal fluctuations in the continuum formulation is necessary to obtain a propercoupling under these demanding constraints.The only other work we are aware of that develops a coupling between a fluctuating compressiblecontinuum solver and a particle method, specifically molecular dynamics, is a coupling schemedeveloped over the last several years by Coveney, Flekkoy, de Fabritiis, Delgado-Buscallioni andcollaborators [2, 43, 44]. There are two important differences between their method and ouralgorithm. Firstly, their method is (primarily, but not entirely) a flux-flux coupling scheme, unlikeour state-flux coupling. Secondly, we do not use MD but rather I-DSMC, which, as discussedbelow, significantly simplifies the handling of the continuum-particle interface.It is not difficult to impose boundary conditions for the continuum subdomain based on theparticle data, and any consistent boundary condition for the PDE being solved can in principlebe imposed. When the continuum solver is deterministic the fluctuations (often inappropriatelyreferred to as “noise”) in the particle data need to be filtered using some sort of spatial and temporalaveraging [41], or the continuum solver needs to be robust [20, 29].The difficult part in coupling schemes is the handling of the boundary conditions for the particlesubdomain. It is very difficult to truncate a fluid particle region without introducing artifacts in thestructure of the particle fluid, and the transport coefficients of the particle fluid (e.g., the pressureand viscosity) critically depend on the fluid structure. In most molecular simulations periodicboundary conditions are used to avoid such artifacts, however, this is not possible in our case.In order to minimize artifacts at the boundary of the particle subdomain, most coupling schemesadd an overlap region in addition to the reservoir region that we described above. In the overlapregion, particles are simulated even though the region belongs to the continuum subdomain. The structure of the fluid in the overlap region is left to adjust to that in the particle subdomain,2thus minimizing the artifacts. At the same time, the dynamics of the particles in the overlapregion is somehow constrained to match the underlying continuum subdomain dynamics using, forexample, various forms of constrained Langevin-type thermostats or (non-holonomic) constrainedMD [42, 45]. To prevent particles from flying outside of the particle subdomain, various artificialconstraining forces are added in the reservoir region, and typically some particle insertion/deletionscheme is used as well. The details can be very different and tricky, and we are not aware of anydetailed comparison or even rudimentary mathematical analysis of the different types of couplingother than the stability analysis of four idealized coupling schemes presented in Ref. [38].We do not use an overlap region in the coupling scheme we present here. Firstly, and mostimportantly, because the I-DSMC fluid is structureless, an overlap region and the associated al-gorithmic complications are simply not required. The addition of an overlap region (added fluidmass) almost always introduces some delay and errors in the coupling, and it is usually assumedthat this surface effect is negligible when the size of the overlap region is small compared to the par-ticle region and when the hydro time steps are much larger than the particle time step. These areassumptions that we do not wish to make as we try to minimize the size of the particle subdomain.The simple and direct coupling scheme we have presented, however, does not work if a structuredstochastic particle fluid is used, such as, for example, a structured stochastic fluid like the StochasticHard-Sphere Dynamics (SHSD) fluid [18]. The SHSD fluid is weakly structured, so that at leastparticle insertion/deletion in the reservoir region is not a problem. However, an overlap region isnecessary to smoothly match the fluid structure at the particle-continuum interface. Constrainingthe dynamics in the overlap region can be done in I-DSMC by introducing additional one-particlecollisions (dissipation) that scatter the particle velocities so that their mean equals the continuumfield values. However, it is not clear how to do this consistently when fluctuations are includedin the continuum solver, and in this paper we restrict our focus on structureless (ideal) stochasticfluids.
III. DETAILS OF THE COUPLING ALGORITHM
The basic ideas behind the state-flux coupling were already described in Section II C. In this sec-tion we describe in detail the two components of the particle-continuum coupling method, namely,the imposition of the continuum state as a boundary condition for the particle subdomain and theimposition of the particle flux as a boundary condition for the continuum subdomain. At the sametime, we will make clear that our coupling is not purely of the state-flux form. The handling of3the continuum subdomain is essentially unchanged from the pure continuum case, with the onlydifference being the inclusion of a refluxing step. The handling of the particle subdomain is morecomplex and explained in greater detail, including pseudocodes for several steps involved in takinga micro time step, including insertion of reservoir particles and the tracking of the particle fluxes.
A. State exchange
We first explain how the state in the reservoir macro cells bordering the particle subdomain,denoted by U ( B ) H , is used by the particle algorithm. The micro cells that are inside the reservoirmacro cells and are sufficiently close to the particle subdomain to affect it during a time-interval∆ t P are labeled as reservoir micro cells . For I-DSMC fluids, assuming the length of the micro cellsalong each dimension is L c (cid:39) D , the micro cells immediately bordering the particle subdomainas well as all of their neighboring micro cells need to be included in the reservoir region. Anillustration of the particle and reservoir regions is given in Figs. 1 and 2.At the beginning of each particle time step reservoir particles are inserted randomly into thereservoir micro cells. The number of particles inserted is based on the target density in the corre-sponding reservoir macro cell. The velocities of the particles are chosen from a Maxwell-Boltzmannor Chapman-Enskog [46] distribution (see discussion in Section III D 2) with mean velocity and tem-perature, and also their gradients if the Chapman-Enskog distribution is used, chosen to match themomentum and energy densities in the corresponding macro cell. The positions of the reservoirparticles are chosen randomly uniformly (i.e., sampled from a Poisson spatial distribution) insidethe reservoir cells which does not introduce any artifacts in the fluid structure next to the interfacebecause the I-DSMC fluid is ideal and structureless. This is a major advantage of ideal fluids overmore realistic structured non-ideal fluids. Note that in Ref. [30] we used a particle reservoir toimplement open boundary conditions in pure particle simulations, the difference here is that thestate in the reservoir cells comes from a continuum solver instead of a pre-specified stationary flowsolution.The reservoir particles are treated just like the rest of the particles for the duration of a particletime step. First they are advectively propagated (streamed) along with all the other particles, andthe total mass, momentum and energy transported by particles advectively through the particle-continuum interface, Φ ( I ) P , is recorded. Particles that, at the end of the time step, are not in eithera reservoir micro cell or in the particle subdomain are discarded, and then stochastic collisions areprocessed between the remaining particles. In the traditional DSMC algorithm, collisions occur4 Figure 2: Illustration of a hybrid simulation of two-dimensional plug flow around a permeable stationarydisk (red). The macroscopic grid is shown (dark-blue lines), with each macro cell composed of 6 × only between particles inside the same micro cell, and thus all of the particles outside the particlesubdomain can be discarded [20, 21]. However, in the I-DSMC algorithm particles in neighboringcells may also collide, and thus the reservoir particles must be kept until the end of the time step.Collisions between a particle in the particle subdomain and a particle in a reservoir cell lead tocollisional exchange of momentum and energy through the interface as well and this contributionmust also be included in Φ ( I ) P . Note that the reservoir particles at the very edge of the reservoirregion do not have an isotropic particle environment and thus there are artifacts in the collisionswhich they experience, however, this does not matter since it is only essential that the particles inthe particle subdomain not feel the presence of the interface.5 B. Flux exchange
After the particle subdomain is advanced for n ex (micro) time steps, the particle flux Φ ( I ) P isimposed as a boundary condition in the continuum solver so that it can complete its (macro)time step. This flux exchange ensures strict conservation, which is essential for long-time stability.Assume that the continuum solver is a one-step explicit method that uses only stencils of widthone, that is, calculating the flux for a given macro cell-cell interface only uses the values in theadjacent cells. Under such a scenario, the continuum solver only needs to calculate fluxes for thecell-cell interfaces between continuum cells, and once the particle flux Φ ( I ) P is known the continuumtime step (5) can be completed. It is obvious that in this simplified scenario the coupling is purelyof the state-flux form. However, in practice we use a method that combines pieces of state and fluxexchange between the particle and continuum regions.First, the particle state is partly used to advance the continuum solver to the next time step.Our continuum solver is a multi-stage method and uses stencils of width two in each stage, thususing an effective stencils that can be significantly larger than one cell wide [9]. While one canimagine modifying the continuum solver to use specialized boundary stencils near the particle-continuum interface (e.g., one-sided differencing or extrapolation), this is not only more complexto implement but it is also less accurate. Instead, the continuum method solves for the fields overthe whole computational domain (continuum patch in AMAR terminology) and uses hydrodynamicvalues for the particle subdomain (particle patch) obtained from the particle solver. These valuesare then used to calculate provisional fluxes Φ ( I ) H and take a provisional time step, as if therewere no particle subdomain. This makes the implementation of the continuum solver essentiallyoblivious to the existence of the particle regions, however, it does require the particle solver toprovide reasonable conserved values for all of the macro particle cells. This may not be possiblefor cells where the continuum hydrodynamic description itself breaks down, for example, cells thatoverlap with or are completely covered by impermeable beads or features of a complex boundary.In such empty cells the best that can be done is to provide reasonable hydrodynamic values, forexample, values based on the steady state compatible with the specified macroscopic boundaryconditions. For partially empty cells, that is, macro cells that are only partially obscured, onecan use the uncovered fraction of the hydro cell to estimate hydrodynamic values for the wholecell. In practice, we have found that as long as empty cells are sufficiently far from the particle-continuum interface (in particular, empty cells must not border the continuum subdomain) theexact improvised hydrodynamic values do not matter much. Note that for permeable beads there6is no problem with empty cells since the fluid covers the whole domain.Secondly, the provisional continuum fluxes are partly used to advance the particle subdomainto the same point in time as the continuum solver. Specifically, a linear interpolation between thecurrent continuum state and the provisional state is used as a boundary condition for the reservoirparticles. This temporal interpolation is expected to improve the temporal accuracy of the coupling,although we are not aware of any detailed analysis. Note that if n ex = 1 this interpolation makesno difference and the provisional continuum fluxes are never actually used.Once the particle solver advances n ex time steps, a particle flux Φ ( I ) P is available and it isimposed in the continuum solver to finalize the provisional time step. Specifically, hydrodynamicvalues in the particle macro cells are overwritten based on the actual particle state, ignoring theprovisional prediction. In order to correct the provisional fluxes, a refluxing procedure is used inwhich the state U ( B ) H in each of the continuum cells that border the particle-continuum interfaceare changed to reflect the particle Φ ( I ) P rather than the provisional flux Φ ( I ) H , U ( B ) H ← U ( B ) H − Φ ( I ) H + Φ ( I ) P . This refluxing step ensures strict conservation and ensures continuity of the fluxes across theinterface, in addition to continuity of the state.
C. Taking a macro time step
Algorithm 1 summarizes the hybrid algorithm and the steps involved in advancing both thesimulation time by one macro time step ∆ t C = n ex ∆ t P . Note that at the beginning of the simula-tion, we initialize the hydrodynamic values for the continuum solver, consistent with the particledata in the particle subdomain and generated randomly from the known (Gaussian) equilibriumdistributions in the continuum subdomain. Algorithm 1:
Take a macro time step by updating the continuum state U H from time t to time t + ∆ t C . Provisionally advance the continuum solver : Compute a provisional macro solution U nextH at time t + ∆ t C everywhere, including the particle subdomain, with an estimated (integrated) provisionalflux Φ ( I ) H through the particle-continuum interface. Reset the particle flux Φ ( I ) P ← .2. Advance the particle solver : Take n ex micro time steps (see Algorithm 2):(a) At the beginning of each particle time step reservoir particles are inserted at the boundary ofthe particle subdomain with positions and velocities based on a linear interpolation between U H and U nextH (see Algorithm 4). This is how the continuum state is imposed as a boundarycondition on the particle subdomain.(b) All particles are propagated advectively by ∆ t P and stochastic collisions are processed, accu-mulating a particle flux Φ ( I ) P (see Algorithm 3).3. Synchronize the continuum and particle solutions : (a) Advance : Accept the provisional macro state, U H ← U nextH .(b) Correct : The continuum solution in the particle subdomain U ( P ) H is replaced with cell averagesof the particle state at time t + ∆ t C , thus forming a composite state over the whole domain.(c) Reflux:
The continuum solution U ( B ) H in the macro cells bordering the particle subdomain iscorrected based on the particle flux, U ( B ) H ← U ( B ) H − Φ ( I ) H + Φ ( I ) P . This effectively imposes theparticle flux as a boundary condition on the continuum and ensures conservation in the hybridupdate.(d) Update the partitioning between particle and continuum cells if necessary. Note that this stepmay convert a continuum cell into an unfilled (devoid of particles) particle cell.
D. Taking a micro time step
Taking a micro time step is described in Algorithm 2, and the remaining subsections give furtherdetails on the two most important procedures used. The initial particle configuration can mosteasily be generated by marking all macro cells in the particle domain as unfilled.
Algorithm 2:
Take an I-DSMC time step. We do not include details about the handling of thenon-DSMC (solute) particles here. Note that the micro time step counter n P should bere-initialized to zero after every n ex time steps.
1. Visit all macro cells that overlap the reservoir region or are unfilled (i.e., recently converted fromcontinuum to particle) one by one, and insert trial particles in each of them based on the continuumstate U H , as described in Algorithm 4.2. Update the clock t ← t + ∆ t P and advance the particle subdomain step counter n P ← n P + 1. Notethat when an event-driven algorithm is used this may involve processing any number of events thatoccur over the time interval ∆ t P [30].3. Move all I-DSMC particles to the present time, updating the total kinetic mass, momentum andenergy flux through the coupling interface F ( I ) P whenever a particle crosses from the particle into thecontinuum subdomain or vice versa, as detailed in Algorithm 3.
4. Perform stochastic collisions between fluid particles, including the particles in the reservoir region.Keep track of the total collisional momentum and energy flux through the coupling interface byaccounting for the amount of momentum ∆ p ij = m ∆ υ ij and kinetic energy ∆ e ij transferred from aregular particle i (in the particle subdomain) to a reservoir particle j (in the continuum subdomain), Φ ( I ) P ← Φ ( I ) P + (cid:2) , ∆ p ij , ∆ e ij (cid:3) .5. Remove all particles from the continuum subdomain.6. Linearly interpolate the continuum state to the present time, U H ← ( n ex − n P ) U H + U nextH n ex − n P + 1
1. Particle flux tracking
As particles are advected during a particle time step, they may cross from a particle to acontinuum cell and vice versa, and we need to keep track of the resulting fluxes. Note that aparticle may cross up to d cell interfaces near corners, where d is the spatial dimension, and mayeven recross the same interface twice near hard-wall boundaries. Therefore, ray tracing is the mostsimple and reliable way to account for all particle fluxes correctly. For the majority of the particlesin the interior, far from corners and hard walls, the usual quick DSMC update will, however, besufficient. Algorithm 3:
Move the fluid particle i from time t to time t + ∆ t P and determine whether itcrosses the coupling interface. We will assume that during a particle time step no particle canmove more than one macro cell length along each dimension (in practice particles typically moveonly a fraction of a micro cell).
1. Store information on the macro cell c old to which the particle belongs, and then tentatively updatethe position of the particle r i ← r i + υ i ∆ t P , and find the tentative macro cell c i and micro cell b i towhich the particle moves, taking into account periodic boundary conditions.2. If c old and c i are near a boundary, go to step 3. If c old ≡ c i , or if c old and c i are both continuum orare both particle cells and at least one of them is not at a corner, accept the new particle positionand go to step 4.3. Undo the tentative particle update, r i ← r i − υ i ∆ t P , and then ray trace the path of the particleduring this time step from one macro cell-cell interface I c to the next, accounting for boundaryconditions (e.g., wrapping around periodic boundaries and colliding the particle with any hard wallsit encounters [30]). Every time the particle crosses from a particle to a continuum cell, update the particle flux at the cell interface I c , Φ ( I c ) P ← Φ ( I c ) P + (cid:2) m, m υ i , mυ i / (cid:3) , similarly, if the particle crossesfrom a continuum to a particle cell, Φ ( I c ) P ← Φ ( I c ) P − (cid:2) m, m υ i , mυ i / (cid:3) .4. If the new particle micro cell b i is neither in the particle subdomain nor in the reservoir region, removethe particle from the system.
2. Inserting reservoir particles
At every particle time step, reservoir particles need to be inserted into the reservoir region orunfilled continuum cells. These particles may later enter the particle subdomain or they may bediscarded, while the trial particles in unfilled cells will be retained unless they leave the particlesubdomain (see Algorithm 3). When inserting particles in an unfilled macro cell, it is importantto maintain strict momentum and energy conservation by ensuring that the inserted particles haveexactly the same total momentum and energy as the previous continuum values. This avoids globaldrifts of momentum and energy, which will be important in several of the examples we present inSection IV. It is not possible to ensure strict mass conservation because of quantization effects, butby using smart (randomized) rounding one can avoid any spurious drifts in the mass.The velocity distribution for the trial particles should, at first approximation, be chosen from aMaxwell-Boltzmann distribution. However, it is well-known from kinetic theory that the presenceof shear and temperature gradients skews the distribution, specifically, to first order in the gradientsthe Chapman-Enskog distribution is obtained [46, 47]. It is important to note that only the kinetic contribution to the viscosity enters in the Chapman-Enskog distribution and not the full viscositywhich also includes a collisional viscosity for dense gases [46]. Previous work on deterministicDSMC hybrids has, as expected, found that using the Chapman-Enskog distribution improvesthe accuracy of the hybrids [20, 41]. However, in the presence of transient fluctuations and theassociated transient gradients it becomes less clear what the appropriate distribution to use is, aswe observe numerically in Section IV A.Certainly at pure equilibrium we know that the Maxwell-Boltzmann distribution is correct de-spite the presence of fluctuating gradients. At the same time, we know that the CE distributionought to be used in the presence of constant macroscopic gradients, as it is required in order toobtain the correct kinetic contribution to the viscous stress tensor [48]. The inability to estimatetime-dependent mean gradients from just a single fluctuating realization forces the use of instanta-neous gradients, which are unreliable due to the statistical uncertainty and can become unphysicallylarge when N is small (say N < a priori (for example, Couette or Fourier flows) and they can be used without trying tonumerically estimate them, otherwise, they can be assumed to be zero or numerically estimated byperforming some spatio-temporal averaging [21]. Note that here we assume that the particles areuniformly distributed (Poisson spatial distribution) as in an ideal gas, and we do not try to takethe density gradient into account (but see the procedure described in the Appendix of Ref. [21]). Algorithm 4:
Insert trial particles in the reservoir or unfilled macro cell c with centroid r c ,taking into account the target density ρ c , velocity v c and temperature T c in cell c , calculatedfrom the conserved cell state U c . If using the Chapman-Enskog distribution, also take intoaccount the estimated local shear rate ∇ c v and local temperature gradient ∇ c T .
1. Build a list L rc of the N rc micro cells contained in the macro cell c that need to be filled with particles.This typically excludes micro cells that are partially covered by impermeable beads or boundaries soas to avoid generating overlaps.2. Determine the total number N p of trial particles to insert into the reservoir portion of c by samplingfrom the binomial distribution P ( N p ) = ¯ N p N p p N p (1 − p ) ¯ N p − N p , where ¯ N p = (cid:98) ρ c V c /m (cid:99) is the total expected number of particles in c , V c is the volume of c , and p = N rc /N sc , where N sc is the number of micro subcells per macro cell. For sufficiently large ¯ N p thiscan be well-approximated by a Gaussian distribution, which can be sampled faster.3. For each of the N p trial particles to be inserted, do:(a) Choose a micro cell b uniformly from the N rc micro cells in the list L rc .(b) Generate a random particle position r i ← r c + r rel uniformly inside micro cell b .(c) Generate a random relative velocity for the particle v rel from the Maxwell-Boltzmann orChapman-Enskog distribution [47], and set the particle velocity by taking into account the de-sired continuum state in cell c and, if available, its estimated gradient, υ i ← v c + ( ∇ c v ) r rel + v rel .(d) If the cell c is an unfilled macro cell, keep track of the total momentum P and energy E of the N p trial particles, to be adjusted in Step 4 for conservation.4. If cell c was unfilled and N p >
1, then correct the particle velocities to match the desired totalmomentum P c = p c V c and energy E c = V c e c inside macro cell c , thus ensuring exact conservation: (a) Calculate the scaling factor α = E c − P c / (2 mN p ) E − P / (2 mN p )and velocity shift ∆ υ = ( P − P c ) / ( αmN p ).(b) Scale and shift the velocity for every trial particle i , υ i ← α ( υ i − ∆ υ ). IV. RESULTS
In this section we provide extensive tests of the hybrid scheme, in both equilibrium and non-equilibrium situations, and in both two and three dimensions. Our goal is to access how well thehybrid method can reproduce results obtained with a pure particle method, which we consider thegold standard.We have implemented our hybrid method in a code that can handle both two and three dimen-sional systems. Of course, one can study one- and two-dimensional flows with the three dimensionalcode by using periodic boundaries along the remaining dimensions. We refer to this as quasi one-or two-dimensional simulations. At the same time, both (I-)DSMC and the LLNS equations havea truly two-dimensional formulation, which we have also implemented for testing purposes. Trans-port coefficients for two-dimensional particle models formally diverge in the infinite-time limit (seethe discussion in Ref. [49]) and it is not obvious that the Navier-Stokes equations are a propercoarse-graining of the microscopic dynamics. However, this divergence is very slow (logarithmic)and it will be mollified (bounded) by finite system size, and we will therefore not need to con-cern ourselves with these issues. In the first three examples, we use the three-dimensional particleand continuum codes, and use the two-dimensional code only for the adiabatic piston example forcomputational reasons.We have also implemented continuum solvers for both the full non-linear and the linearized
LLNS equations. As discussed in more detail in Ref. [9], the nonlinear LLNS equations aremathematically ill-defined and this can lead to breakdown in the numerical solution such as negativedensities or temperatures. At the same time, the linearized equations are not able to describe awide range of physical phenomena such as the effect of fluctuations on the mean flow; they also omita number of terms of order N − , where N is the average number of particles per continuum cell(e.g., the center-of-mass kinetic energy). If the number of particles per continuum cell is sufficientlylarge (in our experience, N >
75) the fluctuations are small and the difference between the linearand nonlinear hydrodynamic solvers is very small, and we prefer to use the nonlinear solver. Wewill use N ∼
100 in our hybrid simulations.2In our implementation, the continuum solver can either be the more accurate third-order Runge-Kutta (RK3) temporal integrator developed in Ref. [9] or a more classical stochastic MacCormackintegrator [6]. The analysis in Ref. [9] shows that obtaining reasonably-accurate equilibriumfluctuations with predictor-corrector methods for the diffusive fluxes, as used in the MacCormackscheme, requires using a continuum time step ∆ t C that is a fraction of the CFL stability limit∆ t CF L . In the simulations we present here we have typically used a time step ∆ t ≈ . t CF L ,which is typically still about 5 times larger than the particle time step ∆ t P , and we have foundlittle impact of the exact value of ∆ t C .The hybrid method requires estimates of the transport coefficients of the particle fluid, notablythe viscosity and the thermal conductivity. For traditional DSMC at low densities there are ratheraccurate theoretical estimates of the viscosity and thermal conductivity [50, 51], however, it isnontrivial to obtain reasonably accurate theoretical values at the higher densities we use in I-DSMC because of the importance of multi-particle correlations. We therefore estimate the transportcoefficients numerically. For this purpose we simulate a system with periodic boundaries along twoof the directions and isothermal stick wall boundaries along the other direction. To estimate theviscosity we apply a shear flow by moving one of the wall boundaries at constant speed inducinga Couette flow with an approximately constant shear gradient ¯ ∇ v = ( ∇ v + ∇ v T ). We thencalculate the steady-state stress tensor σ by averaging over all particles i and colliding pairs ofparticles ij over a long time interval ∆ t , σ = σ k + σ c = m (cid:104) υ i ⊗ υ i (cid:105) + (cid:10) r ij ⊗ ∆ p ij (cid:11) c ∆ t , where σ k = P I +2 η k ¯ ∇ v is the kinetic contribution giving the kinetic viscosity η k , and σ c = 2 η c ¯ ∇ v is the collisional contribution giving the collisional viscosity η c , η = η k + η c . We exclude particlesthat are close to the wall boundaries when calculating these averages to minimize finite-size effects.Similarly, for the thermal conductivity we apply a small constant temperature gradient ∇ T bysetting the two walls at different temperatures, and we also impose the required density gradientto maintain mechanical equilibrium (constant pressure). We then calculate the steady-state heatflux vector ξ = µ ∇ T , ξ = ξ k + ξ c = m (cid:28) υ i υ i (cid:29) + (cid:104) (∆ e ij ) r ij (cid:105) c ∆ t , from which we obtain the kinetic and collisional contributions to the thermal conductivity. Thereare alternative methods that one can use to calculate the transport coefficients, using both equilib-rium and non-equilibrium settings, however, we have found the above method to be most accurate3for a given computational effort if only moderate accuracy is desired. Results from different non-equilibrium methods are found to be within 5 −
10% of each other.
A. Mismatch at the interface
Previous work has studied a hybrid scheme very similar to the one described here for severalquasi one-dimensional situations [21]. Reference [21] first studied a pure equilibrium situation inwhich one part of a periodic domain was covered by a particle subdomain, and found that thestochastic hybrid scheme was able to reproduce the spatio-temporal correlations in equilibriumfluctuations very well, although some mismatch at the particle/continuum interface was found.Here we explore this mismatch more carefully by studying a quasi one-dimensional periodic systemwhere the middle portion of the domain is filled with particles and the rest is continuum. In eachmacro cell, we compute the average density, temperature and velocity and their variances withhigh accuracy.We first calculate the mean conserved quantities (mean density, momentum density, and energydensity) in each macro cell, and then calculate the mean velocity and temperature from those, forexample, (cid:104) v (cid:105) = (cid:104) j (cid:105) / (cid:104) ρ (cid:105) instead of averaging the instantaneous velocities, (cid:104) v (cid:105) = (cid:104) j/ρ (cid:105) . As shownin Ref. [52], the approach we use leads to an unbiased estimate of the mean, while the latter has abias when there are correlations between the fluctuations of the different hydrodynamic variables.The variances are estimated from the instantaneous values, e.g., (cid:10) δv (cid:11) = (cid:68) ( j/ρ ) (cid:69) − (cid:104) j/ρ (cid:105) . It ispossible to construct unbiased estimates for the variances as well [52], however, this is somewhatmore involved and the bias is rather small compared to the artifacts we are focusing on.In Fig. 3 we show the means and variances along the length of the system, normalized bythe expected values [6]. For velocity, the mean is zero to within statistical uncertainty in boththe particle and continuum subdomains, and we do not show it in the figure. However, a smallmismatch is clearly seen in Fig. 3(a) between the density and temperature in the particle andcontinuum subdomains. The mismatch is such that the pressure (cid:104) p (cid:105) = (cid:104) ρ (cid:105) R (cid:104) T (cid:105) is constant acrossthe interface, that is, the particle and continuum subdomains are in mechanical equilibrium but notin true thermodynamic equilibrium. In Appendix A we show that this kind of mismatch is expectedbecause the average particle fluxes coming from the reservoir particles inserted in the continuumsubdomains have a bias of order N − , where N is the average number of particles per macroscopiccell. Because our coupling matches both the state and the fluxes across the interface this biasmakes it impossible for the particle and continuum to reach true thermodynamic equilibrium. The4theory in Appendix A suggests that the size of the mismatch is of order N − , consistent withthe results in Fig. 3(a); however, the crude theoretical estimates do not actually give the steadystate since they assume equilibrium to begin with. That the cells near the interface are not inequilibrium is reflected in the variances of the hydrodynamic variables, which have a spike near theparticle-continuum interface, as shown in Fig. 3(b). We have observed that the relative magnitudeof this spike does not depend on N .The cause of the mismatch is the fact that we use the instantaneous values of the local density,velocity and temperature when generating the velocities for the reservoir particles. This is necessarybecause we cannot in general obtain estimates of the time-dependent mean values from running asingle realization, and are forced to use the instantaneous values. One can use some sort of spatio-temporal averaging to obtain estimates of the local means, however, this introduces additionaltime and length scales into the algorithm that do not have an obvious physical interpretation. Forsteady-state problems it may be possible to use running means to avoid the mismatch we observe,however, this is not possible in a general dynamic context. Deeper theoretical understanding ofthe connection between the microscopic dynamics and the LLNS equations is necessary to designa more consistent approach.Another complex issue that arises in the fluctuating hybrid is whether to use the Maxwell-Boltzmann (MB) or the Chapman-Enskog (CE) distributions when generating the velocities ofthe reservoir particles (see discussion in Section III D 2). In Fig. 3(a) we compare the size of themismatch at the interface when the MB and CE distributions are used. For the CE distributionwe obtained a local estimate of the gradient using simple centered differences of the instantaneoushydrodynamic variables. As seen in the figure, there is a greater discrepancy when the CE dis-tribution is used, especially for the variances. We will therefore adopt a compromise in which weuse the MB distribution for all calculations reported here. In cases when there is a macroscopicbackground gradient that is specified a priori (e.g., shear flows) we can use that gradient in theCE distribution.Note that we have performed numerous additional quasi-one dimensional tests of the couplingthat we do not describe here for brevity. For example, we have tested the matching of the shearstress tensor in a shear flow parallel to the particle-continuum interface by verifying that a linearvelocity profile is obtained without any slope discontinuity at the interface (see Appendix D.2 inRef. [44]). We have also reproduced the results reported in Ref. [21], such as the presence ofunphysical long-range correlations in the fluctuations in the particle region when the deterministicinstead of the stochastic hybrid is used.5 x N o r m a li z e d m e a n ρ (MB, N=120)T (MB, N=120) ρ (CE, N=120)T (CE, N=120) ρ (MB, N=480) T (MB, N=480) (a) x N o r m a li z e d v a r i a n ce v x (MB)v y (MB) ρ (MB)T (MB)v x (CE)v y (CE) ρ (CE)T (MB) (b) Figure 3: Normalized means and variances of the hydrodynamic variables in each of the 46 macroscopiccells of a quasi one-dimensional periodic system where the middle portion (20 continuum cells, in-betweendashed vertical lines) is the particle subdomain. The velocities of the reservoir particles are either samplesfrom a Maxwell-Boltzmann (MB, lines only) or a Chapman-Enskog distribution (CE, lines and symbols).(a) Unbiased [52] estimates of the mean density and temperature (we averaged about 2 . · samples takenevery macro step) for two different sizes of the continuum cells, ones containing N = 120 particles onaverage, and ones containing N = 480 particles. Only the MB results are shown for the larger cells forclarity, with similar results observed for CE. (b) Estimates of the variances of the hydrodynamic variablesfor N = 120 particles per continuum cell (also an average over 2 . · samples). B. Dynamic Structure Factors
The hydrodynamics of the spontaneous thermal fluctuations in the I-DSMC fluid is expectedto be described by the Landau-Lifshitz Navier-Stokes (LLNS) equations for the fluctuating field U = ( ρ + δρ, δv, T + δT ) linearized around a reference equilibrium state U = ( ρ , v = , T )[9, 35, 53]. By solving these equations in the Fourier wavevector-frequency ( k , ω ) domain for (cid:98) U ( k , ω ) = ( (cid:98) δρ, (cid:99) δv, (cid:99) δT ) and performing an ensemble average over the fluctuating stresses we canobtain the equilibrium (stationary) spatio-temporal correlations (covariance) of the fluctuatingfields. We express these correlations in terms of the 3 × hydrodynamicstructure factor matrix S H ( k , ω ) = (cid:68) (cid:98) U (cid:98) U (cid:63) (cid:69) [9], which is essentially the spatio-temporal spectrumof the fluctuating fields. Integrating S H ( k , ω ) over frequencies gives the hydrostatic structure factormatrix S H ( k ), which turns out to be diagonal since in any given snapshot the hydrodynamicvariables are uncorrelated at equilibrium.We non-dimensionalize S H ( k , ω ) so that S H ( k ) is the identity matrix. For example, the density-density correlations are given by the dimensionless structure factor S ρ ( k , ω ) = (cid:0) ρ c − k B T (cid:1) − (cid:104) ˆ ρ ( k , ω ) ˆ ρ (cid:63) ( k , ω ) (cid:105) , and we express the spatio-temporal cross-correlation between density and velocity through thedimensionless structure factor S ρ,v ( k , ω ) = (cid:0) ρ c − k B T (cid:1) − (cid:0) ρ − k B T (cid:1) − (cid:104) ˆ ρ ( k , ω )ˆ v (cid:63) ( k , ω ) (cid:105) , where c = k B T /m is the isothermal speed of sound. Reference [21] demonstrated that a hybridmethod very similar to the one we described here correctly reproduces the density-density time-correlation function S ρ ( k , t ) for large wavelengths in a quasi one-dimensional periodic system. Thedensity-density dynamic structure factor S ρ ( k , ω ) is often the only one considered because it isaccessible experimentally via light scattering measurements and thus most familiar [53]. The fulldynamic structure factor matrix S H ( k , ω ) is a more complete measure of the spatio-temporal evo-lution of the thermal fluctuations and includes both sound (hyperbolic) and dissipative (diffusive)effects. It is therefore important to show that the hybrid scheme correctly reproduces S H ( k , ω ) ascompared to a purely particle simulation and demonstrate that the hybrid is capable of capturingthe propagation of spontaneous thermal fluctuations across the particle-continuum interface.7
1. Bulk Dynamic Structure Factor
For a bulk fluid, i.e., for periodic boundary conditions, it is well-known [53] that the density-density component S ρ ( k , ω ) and the temperature-temperature component S T ( k , ω ) of S H ( k , ω )exhibit three peaks for a given wavevector k . There is one central Rayleigh peak at ω = 0 whichcomes from entropic fluctuations. There is also two symmetric Brillouin peaks at ω ≈ c s k , where c s is the adiabatic speed of sound, which come from the isoentropic propagation of sound wavesinduced by the fluctuations. For the components of the velocity parallel to the wavevector thedynamic structure factors S v (cid:107) ( k , ω ) exhibit all three peaks, while for the component perpendicularto the wavevector S v ⊥ ( k , ω ) lacks the central peak.Figure 4 shows selected dynamic structure factors for a quasi two-dimensional system withperiodic boundary conditions. The simulation box is composed of 10 × × x axes and designatedas a particle subdomain, and the remaining two thirds of the domain were continuum.In Fig. 4 we show the results for a wavevector k that is neither parallel nor perpendicular tothe particle-continuum interface so as to test the propagation of both perpendicular and tangen-tial fluctuations across the interface. The results show very little discrepancy between the pureparticle and the hybrid runs, and they also conform to the theoretical predictions based on theLLNS equations. Perfect agreement is not expected because the theory is for the spectrum of thecontinuum field while the numerical results are discrete spectra of cell averages of the field, a dis-tinction that becomes important when the wavelength is comparable to the cell size. Additionally,even purely continuum calculations do not reproduce the theory exactly because of spatio-temporaldiscretization artifacts.
2. Dynamic Structure Factors for Finite Systems
The previous section discussed the “bulk” dynamic structure factors, as obtained by using peri-odic boundary conditions. For non-periodic (i.e., finite) systems equilibrium statistical mechanics8 -20 -15 -10 -5 0 5 10 15 20 ω S ( k , ω ) ρ Hybridv x v y T ρ Theory ρ Particle -20 -15 -10 -5 0 5 10 15 20 ω -0.5-0.4-0.3-0.2-0.100.10.2 S ( k , ω ) ρ − v x Hybridv x - v y ρ − Tv x - T ρ − v x Theory ρ − v x Particle
Figure 4: Discrete dynamic structure factors for waveindices q x = k x L x / (2 π ) = 1, q y = 2 and q z = 0 for aquasi two-dimensional system with lengths L x = L y = 2, L z = 0 .
2, split into a grid of 10 × × D = 0 .
04 (density φ = 0 .
5, and collisionfrequency χ = 0 . temporal snapshots is employed. Weperform pure particle runs and also hybrid runs in which only a strip 4 macro cells along the x axes was filledwith particles and the rest handled with a continuum solver. The different hydrodynamic pairs of variablesare shown with different colors, using a solid line for the result from the hybrid runs, symbols for the resultsof the pure particle runs, and a dashed line for the theoretical predictions based on the linearized LLNSequations (solved using the computer algebra system Maple). ( Top ) The diagonal components S ρ ( k , ω ), S v x ( k , ω ), S v y ( k , ω ) and S T ( k , ω ). ( Bottom ) The off-diagonal components (cross-correlations) S ρ,v x ( k , ω ), S v x ,v y ( k , ω ), S ρ,T ( k , ω ) and S v x ,T ( k , ω ). ∂ Ω with normal vector n either Dirichlet or von Neumannboundary conditions need to be imposed on the components of the velocity and the temperature(the boundary condition for density follows from these two). Two particularly common types ofboundaries are: Thermal walls for which a stick condition is imposed on the velocity, v ∂ Ω = 0, and the temper-ature is fixed, T ∂ Ω = T . Adiabatic walls for which a slip condition is imposed on the velocity, n · ∇ v (cid:107) = 0 and v ⊥ = 0,and there is no heat conduction through the wall, n · ∇ T = 0.In particle simulations, these boundary conditions are imposed by employing standard rules forparticle reflection at the boundaries [30]. We describe the corresponding handling in continuumsimulations in Appendix C. In Appendix B we derive the form of the additional peaks in thedynamic structure factor for adiabatic walls by solving the linearized LLNS equations with theappropriate conditions.In Figure 5 we show dynamic structure factors for a quasi one-dimensional system bounded byadiabatic walls. As for the bulk (periodic) case in the previous section, we perform both purelyparticle runs and also hybrid runs in which the middle third of the domain is designated as aparticle subdomain. Additional peaks due to the reflections of sound waves from the boundariesare clearly visible and correctly predicted by the LLNS equations and also accurately reproducedby both the purely continuum solution (not shown) and the hybrid. Similar agreement (not shown)is obtained between the particle, continuum and the hybrid runs for thermal walls. These resultsshow that the hybrid is capable of capturing the dynamics of the fluctuations even in the presenceof boundaries. Note that when the deterministic hybrid scheme is used one obtains essentially thecorrect shape of the peaks in the structure factor (not shown), however, the magnitude is smaller(by a factor of 2 . C. Bead VACF
As an illustration of the correct hydrodynamic behavior of the hybrid algorithm, we studythe velocity autocorrelation function (VACF) C ( t ) = (cid:104) v x (0) v x ( t ) (cid:105) for a large neutrally-buoyant ω S ( k , ω ) Theory (adiabatic)ParticleHybridPeriodic S(k, ω )Periodic S(k, ω )/2 ω S ( k , ω ) Theory (adiabatic)ParticleHybridPeriodic S(k, ω )Periodic S(k, ω )/2 ω Figure 5: Discrete dynamic structure factors for a quasi one-dimensional system with length L = 7 . q = 2 and wavevector k = 2 πq/L ≈ .
75. The I-DSMC fluid parameters are as in Fig. 4. We perform pure particle runs andalso hybrid runs in which the middle third of the domain is filled with particles and the rest handled witha continuum solver. The results from purely particle runs are shown with symbols, while the results fromthe hybrid are shown with a solid line. The theoretical predictions based on the linearized LLNS equationsand the equations in Appendix B (solved using the computer algebra system Maple) are shown with adashed line for adiabatic boundaries and dotted line for periodic boundaries. Since the magnitude of theBrillouin peaks shrinks to one half the bulk value in the presence of adiabatic walls, we also show the resultwith periodic boundaries scaled by 1 / Top ) Dynamic structure factor for density, S ρ ( k, ω ), showing the Rayleigh peak and the multiple Brillouin peaks. ( Bottom ) Dynamic structure factorfor the component of velocity perpendicular to the wall, S v ⊥ ( k, ω ), which lacks the Rayleigh peak. Thecorresponding correlations for either of the parallel velocity components, S v (cid:107) ( k, ω ), have only a Rayleighpeak, shown in the inset on a log-log scale. bead of mass M and radius R diffusing through a dense Maxwell I-DSMC stochasticfluid [18] of particles with mass m (cid:28) M and collision diameter D (cid:28) R and density (volumefraction) φ [mass density ρ = 6 mφ/ ( πD )]. The VACF problem is relevant to the modelingof polymer chains or (nano)colloids in solution (i.e., complex fluids ), in particular, the integralof the VACF determines the diffusion coefficient which is an important macroscopic quantity.Furthermore, the very first MD studies of the VACF for fluid molecules led to the discovery of along power-law tail in C ( t ) [54] which has since become a standard test for hydrodynamic behaviorof methods for complex fluids [55, 56, 57, 58, 59, 60, 61].The fluctuation-dissipation principle [62] points out that C ( t ) is exactly the decaying speed of abead that initially has a unit speed, if only viscous dissipation was present without fluctuations, andthe equipartition principle tells us that C (0) = (cid:10) v x (cid:11) = kT / M . Using these two observations andassuming that the dissipation is well-described by a continuum approximation with stick boundaryconditions on a sphere of radius R H , C ( t ) has been calculated from the linearized (compressible)Navier-Stokes (NS) equations [63, 64]. The results are analytically complex even in the Laplacedomain, however, at short times an inviscid compressible approximation applies. At large timesthe compressibility does not play a role and the incompressible NS equations can be used to predictthe long-time tail [64, 65]. At short times, t < t c = 2 R H /c s , the major effect of compressibilityis that sound waves generated by the motion of the suspended particle carry away a fraction ofthe momentum with the sound speed c s , so that the VACF quickly decays from its initial value C (0) = k B T /M to C ( t c ) ≈ k B T /M eff , where M eff = M +2 πR ρ/ t > t visc =4 ρR H / η , the VACF decays as with an asymptotic power-law tail ( k B T /M )(8 √ π ) − ( t/t visc ) − / ,in disagreement with predictions based on the Langevin equation (Brownian dynamics), C ( t ) =( k B T /M ) exp ( − πR H ηt/M ).We performed purely particle simulations of a diffusing bead in various I-DSMC fluids in Refs.[18, 61]. In purely particle methods the length of the runs necessary to achieve sufficient accuracyin the region of the hydrodynamic tail is often prohibitively large for beads much larger than thefluid particles themselves. It is necessary to use hybrid methods and limit the particle region to thevicinity of the bead in order to achieve a sufficient separation of the molecular, sonic, viscous, anddiffusive time scales and study sufficiently large box sizes over sufficiently long times. In the resultswe report here we have used an impermeable hard bead for easier comparison with existing theory;similar results are obtained using permeable beads. The interaction between the I-DSMC fluidparticles and the bead is treated as if the fluid particles are hard spheres of diameter D s , chosento be somewhat smaller than their interaction diameter with other fluid particles (specifically,2we use D s = D/
4) for computational efficiency reasons, using an event-driven algorithm [30].Upon collision with the bead the relative velocity of the fluid particle is reversed in order toprovide a no-slip condition at the surface of the suspended sphere [30, 58] (slip boundaries givequalitatively identical results). We have estimated the effective (hydrodynamic) colloid radius R H from numerical measurements of the Stokes friction force F = − πR H ηv and found it to besomewhat larger than the hard-core collision radius R + D s /
2, but for the calculations below weuse R H = R + D s /
2. Since we used periodic boundary conditions with a box of length L , there areartifacts in C ( t ) after about the time at which sound waves generated by its periodic images reachthe particle, t L = L/c s . The averaging procedure that we used in order to eliminate some of thenoise from the tail does not properly resolve these sound effects , even though they are visible inthe results shown below. Furthermore, there are strong finite-size effects that manifest themselvesas a rapid decay of the VACF for times longer than the viscous time-scale ρL /η [59].For the hybrid calculations, we localize the particle subdomain to the continuum cells thatoverlap or are close to the moving bead. The location of the particle subdomain is updatedperiodically as the bead moves. The algorithm that we use tries to fit the particle subdomain asclosely around the bead but ensuring that there are a certain number of micro cells in-between thesurface of the bead and the particle-continuum interface. The exact shape of the particle subdomainthus changes as the bead moves and the number of particles employed by the hybrid fluctuates,especially when the bead is small compared to the continuum cells. Obtaining reasonably-accurateresults for the VACF at long times requires very long runs. We found that it is crucial to strictlyconserve momentum in the hybrid when unfilled continuum cells are transformed into particlecells. Otherwise very slow drifts in the momentum of the system appear due to the use of periodicboundary conditions, and this drift changes the tail of C ( t ), especially for massive beads wherethe typical bead speed is already small compared to the typical fluid particle speed. We used theMacCormack solver and the linearized formulation of the LLNS equations for these simulations,however, similar results are obtained with the non-linear Runge-Kutta solver as long as the macrocells are sufficiently large. As discussed earlier, the use of the linearized formulation makes itpossible to reduce the size of the continuum cells without introducing numerical problems due to To calculate the VACF, we first calculated the average mean-square displacement ∆ r ( t ) by averaging over a longrun and numerically differentiated this to obtain a time-dependent diffusion coefficient D ( t ). We then smoothed D ( t ) by fitting a quadratic polynomial over short-time intervals spaced on a logarithmic scale (so that more pointswere averaged over in the tail), and obtained the VACF by differentiating the smoothed D ( t ). This procedureproduces well-resolved tails, however, it obscures features over short time scales compared to the smoothing interval.A more direct velocity autocorrelation calculation was used at very short times. t / t visc M C ( t ) / k B T Stoch. hybrid (L=1)Det. hybrid (L=1)Stoch. hybrid (L=2)Det. hybrid (L=2)Particle (L=1)Theory t c s / R t L=1 (a) t / t visc M C ( t ) / k B T Stoch. hybrid (L=2)Det. hybrid (L=2)Stoch. hybrid (L=3)Det. hybrid (L=3)Particle (L=2)Theory t c s / R t L=2 (b)
Figure 6: Normalized velocity autocorrelation function (VACF) C ( t ) / ( k B T /M ) for a neutrally-buoyant hardspherical bead of mass M suspended in a fluid of I-DSMC particles of diameter D , for two different beadsizes, a small bead of radius R = 1 . D (a) and a large bead of radius R = 6 . D (b). A log-log scale is usedto emphasize the long-time power law tail and the time is normalized by the viscous time t visc , so that theresults should be approximately independent of the actual bead radius. The inset shows the initial decay ofthe VACF on a semi-log scale, where the time is normalized by the sonic time scale t c . Periodic boundaryconditions with a cubic cell of length L are employed, and the sound crossing time t L is indicated. Resultsfrom purely particle simulations are shown with a dashed-dotted line, and the incompressible hydrodynamictheory is shown in with a dotted line. Results from hybrid runs are also shown with a solid line for thestochastic hybrid and dashed line for the deterministic hybrid, for both the same box size as the particleruns (red) and a larger simulation box (green). In the calculations reported in Fig. 6, the I-DSMC fluid has a density (volume fraction) φ = 0 . χ = 0 .
62. The adiabatic sound speed is c s = (cid:112) k B T / (3 m ) andviscosity is η = ˜ ηD − √ mk B T , where we measured ˜ η ≈ .
75. Note that in atomistic time units t = D (cid:112) m/k B T the viscous time scale is t visc /t ≈ φ ( R H /D ) / (3 π ˜ η ) ≈ . R H /D ) .As a first test case, in Fig. 6(a) we compare against the particle data from Ref. [61], for whichthe size of the bead is R = 1 . D , M = 7 . m , and the simulation box is L = 1 = 25 D , whichcorresponds to 24 micro cells and about N ≈ . · particles. The hybrid runs used macrocells each composed of 4 micro cells, which corresponds to about N = 80 particles per cell, andthe size of the particle subdomain fluctuated between about 3 · and 6 · particles due to thechange of the location of the bead relative to the continuum grid. The particle result in Fig. 6(a)is the average over 10 runs, each of length T /t visc ≈ · , while the hybrid results are from a4single run of length T /t visc ≈ . · . It is seen in the figure that both the deterministic and thefluctuating hybrid reproduce the particle results closely, with a small but visible difference at longtimes where the deterministic hybrid under-predicts the magnitude of the tail in the VACF. Wealso show results from a hybrid run with a twice larger simulation box, L = 2, which marginallyincreases the computational effort in the hybrid runs, but would increase the length of purelyparticle runs by an order of magnitude. The hydrodynamic tail becomes pronounced and closer tothe theoretical prediction for an infinite system, as expected. We have observed (not shown) thatusing small continuum cells composed of only 3 micro cells, which corresponds to about N = 35particles per cell, leads to an over-prediction of the magnitude of the VACF at short times, thatis, to an excess kinetic energy for the bead by as much as 20%, depending on the exact parametersused.As a second, more difficult test, in Fig. 6(b) we report results from particle simulations fora much larger bead, R = 6 . D , M = 976 m , and the simulation box is L = 2 = 50 D , whichcorresponds to N ≈ . · particles. We have performed a variety of hybrid runs and in Fig.6(b) we report results from runs with macro cells each 3 micro cells, as well as results for a largersimulation box, L = 3, and macro cells each composed of 4 micro cells. The particle results arethe average over 5 runs each of length T /t visc ≈ . · , while the hybrid results are from a singlerun of length T /t visc ≈ . · . The hybrid runs had a particle subdomain containing about 2 · particles. We observed little impact of the size of the continuum cell size or the size of the particlesubdomain. The results in Fig. 6(b) show that the stochastic hybrid correctly reproduces the tailin the VACF, while it slightly over-estimates the VACF at short times. The deterministic hybrid,on the other hand, strongly under-estimates the magnitude of C ( t ) at both short and long times.It is particularly striking that the deterministic hybrid fails to reproduce the magnitude of the longtime tail (and thus the diffusion coefficient), vividly demonstrating the importance of includingfluctuations in the continuum domain.Computing the VACF for a diffusing bead has become a standard test for micro- and nano-scalefluid-structure coupling methods and has been performed for a suspended bead in a wide range ofparticle and (semi-)continuum compressible and incompressible fluids [54, 55, 56, 57, 58, 59, 60, 61].However, these tests often do not correctly test for all of the components necessary to match theVACF over all relevant timescales: equipartition, a fluctuation-dissipation relation, and hydrody-namics. Purely continuum fluid methods allow for using a much larger time step than particle(and thus hybrid) methods, especially if an incompressible formulation is directly coupled to theequations of motion of the suspended bead [57, 59, 60]. When fluctuations are not included in5continuum methods, the VACF is often obtained by considering the deterministic decay of the ve-locity of a bead. This however assumes a priori that proper thermodynamic equilibrium exists withthe correct fluctuation-dissipation relation, without actually testing this explicitly. Alternatively,coupling methods between a continuum fluid and a suspended particle often have some arbitrarycoupling parameters that are tuned to reproduce the desired diffusion coefficient without producinga physically-consistent VACF, especially at short times [66]. In particular, incompressible formu-lations cannot reproduce the initial value or the decay of the VACF and should instead aim toproduce an average kinetic energy of the bead of k B T rather than the statistical mechanics resultof 3 k B T /
D. Adiabatic Piston
The problem of how thermodynamic equilibrium is established has a long history in statisticalmechanics [67]. The adiabatic piston problem is one of the examples used to study the fluctuation-driven relaxation toward equilibrium [68, 69, 70, 71] that is simple enough to be amenable totheoretical analysis but also sufficiently realistic to be relevant to important problems in nano-science such as Brownian motors [24, 72]. We study the following formulation of the adiabaticpiston problem. A long quasi one-dimensional box with adiabatic walls is divided in two halves witha thermally-insulating piston of mass M (cid:29) m that can move along the length of the box withoutfriction. Each of the two halves of the box is filled with a fluid that is, initially, at a differenttemperature T and density ρ , here assumed to follow the ideal equation of state P = ρk B T /m .If the macroscopic pressures in the two halves are different, ρ L T L (cid:54) = ρ R T R , then the pressuredifference will push the piston to perform macroscopic oscillations with a period that can beestimated by assuming that each half undergoes an adiabatic transformation ( P V γ = const.).These oscillations are damped by viscous friction and lead to the piston essentially coming to rest ina state of mechanical equilibrium , ρ L T L = ρ R T R . This stage of the relaxation from non-equilibriumto mechanical equilibrium has been shown to be well-described by deterministic hydrodynamics[70].The state of mechanical equilibrium is however not a state of true thermodynamic equilibrium,which also requires equality of the temperatures on the two sides of the piston. Reaching fullequilibrium requires heat transfer through the piston, but the piston is adiabatic and does not6conduct heat. In classical deterministic hydrodynamics the piston would just stand still and neverreach full equilibrium. It has been realized long ago that heat is slowly transferred through themechanical asymmetric fluctuations of the piston due to its thermal motion, until the temperatureson both sides of the piston equilibrate and the fluctuations become symmetric. This equilibrationtakes place through a single degree of freedom (the piston position) coupling the two large reservoirs,and it would be astronomically slow in a laboratory setting. While various Langevin or kinetictheories have been developed for the effective heat conduction of the adiabatic piston (see Refs.[68, 69, 70, 71] and references therein), there is no complete theoretical understanding of theeffective heat conductivity, especially in dense fluids. Molecular dynamic simulations have beenperformed in the past [68, 69, 70] using hard-disk fluids, but the very long runs required to reachthermodynamic equilibrium for massive pistons have limited the size of the systems that couldbe studied. Furthermore, there have been no studies applying fluctuating hydrodynamics to thisproblem.Here we apply our hybrid method to the adiabatic piston problem in two dimensions, using anon-linear two-dimensional implementation of the Runge-Kutta integrator described in Ref. [9] asthe continuum solver. The choice of two dimensions is for purely computational reasons. Firstly,the number of particles required to fill a box of sufficient size is much smaller thus allowing for longparticle simulations. Secondly, in order to implement the piston in our particle scheme we reusedthe same mixed event-driven/time-driven handling [30] as we used for the VACF computations inSection IV C. Namely, we made a piston out of N b small impermeable beads, connected togetherto form a barrier between the two box halves, as illustrated in Fig. 7. In two dimensions, byensuring that two piston beads never separate by more than a given distance we can ensure thattwo I-DSMC particles on opposite sides of the piston cannot possibly collide and thus the pistonwill be insulating. We have studied two different types of pistons, a flexible piston where thethe beads are tethered together to form a chain [30] that is stretched but where each individualbead can still move independently of the others, and a rigid piston that is obtained with a slightmodification of the event loop to move all of the piston beads in unison. While at the macroscopiclevel the exact shape of the piston should not make a big difference, we have found that increasingthe number of degrees of freedom of the piston from one to N b makes a significant differencein the thermal conductivity of the piston, and therefore, we will focus here on rigid pistons asin the traditional formulation. We use specular collisions of the fluid particles with the pistonbeads, although qualitatively identical (but not quantitatively identical) results are obtained usingbounce-back collisions as well.7 Figure 7: An illustration of the computational setup used for the adiabatic piston computations. Only thecentral portion of the box of aspect ratio 6 × The hybrid method setup for the adiabatic piston is illustrated in Fig. 7. We use a two-dimensional Maxwell I-DSMC particle fluid ( φ = 1, χ = 1) with collision diameter D = 0 . D s = D/
2) and a piston composed of N b = 40 beads of diameter D b = 0 . y dimension (parallel to the piston) with the width ofthe domain L y = 4 being 40 microscopic cells, while adiabatic walls were placed at the ends ofthe box whose total length L x = 24 was 240 microscopic cells. We have studied various sizes forthe macroscopic cells, and report results for a quasi one-dimensional continuum grid in which eachmacro cell contains 4 ×
40 micro cells, corresponding to about 200 particles per continuum cell. Wealso present results for a two-dimensional continuum grid where each macro cell contains 8 × M = 4000 m that started at a position x = 8 that is8 t x ( t ) ParticleStoch. hybridDet. hybrid0 250 500 750 100066.577.58
Figure 8: Relaxation of a massive rigid piston (
M/m = 4000) from a position x = 8 that is not in mechanicalequilibrium. Through rapid damped oscillations mechanical equilibrium is established at position x ≈ not in mechanical equilibrium and thus performs rapid damped oscillations until it reaches a state ofmechanical equilibrium at x ≈
7, from which it slowly relaxes toward true equilibrium. The resultsin the figure show that the stochastic hybrid reproduces the correct relaxation toward equilibriumwhile the deterministic hybrid severely under-predicted the rate of equilibration (effective heatconductivity), even though the initial mechanical stage of the relaxation is correctly captured byboth hybrids, as expected. We have observed that the deterministic hybrid fails to give the correctanswer whenever a rigid massive piston is used,
M > m . For flexible pistons, we find that evenfor a large bead mass M b (overall piston mass M = N b M b ) both the deterministic and stochastichybrids reproduce the purely particle results for the slow relaxation toward equilibrium.A more detailed comparison of the particle and hybrid results for a piston of mass M = 1000 m that is initially in mechanical equilibrium at position x = 8 = L x / k B T L = 2 / ρ L = 2 / k B T R = 4 / ρ L = 1 /
3, so that there is an equal masson each side of the piston. At the true equilibrium state the piston remains close to the middle9 × × × × t x ( t ) Stoch. hybrid (4x40)Det. hybrid (4x40)Stoch. hybrid (8x10)Particle0 5 × × P i s t o n o ff s e t Figure 9: Relaxation of a rigid piston of mass
M/m = 1000 from an initial state of mechanical equilibrium( x = 8) to a state of thermodynamic equilibrium ( x = 12). The inset emphasizes the initial exponentialdecay on a semi-log scale. The hybrid runs used a particle subdomain of width w P = 2 on each side ofthe piston and continuum cells that were composed of either 4 ×
40 or 8 ×
10 microscopic cells. For thedeterministic hybrid the macro cell size makes little difference so we only show the 4 ×
40 case. of the box, x eq = L/ position that are diminished by direct averaging because they have different (random) phases.Figure 9 shows that the stochastic hybrid is able to correctly reproduce the rate of exponentialrelaxation of the piston toward equilibrium with many fewer particles than the purely particle runs,while the deterministic hybrid fails. We have observed a slight dependence on the exact details ofthe hybrid calculations such as cell size or the width of the particle subdomain; however, in general,the stochastic hybrid has shown to be remarkably robust and successful. At the same time, theimportance of including thermal fluctuations in the continuum subdomain is revealed as for theVACF computations in Section IV C.The relation to the equilibrium VACF computations in Section IV C is emphasized by computingthe VACF C ( t ) = (cid:104) v x ( t ) v x (0) (cid:105) for the piston in its state of true equilibrium , k B T = 1. The At thermodynamic equilibrium with a common temperature T , the frequency of the thermally-driven oscillationscan be estimated using a quasi-adiabatic harmonic approximation to be ω ≈ Nk B T / [ Mx ( L x − x )] (see also analternative derivation in Ref. [69]) and the amplitude of the oscillations can be estimated to be on the order of∆ x ≈ x ( L x − x ) /N , where N is the number of fluid particles per chamber. In Ref. Ref. [69] the autocorrelation function for the piston bead position is used to extract the rate of exponential × -4 × -3 ParticleStoch. w P =2 Det. w P =2 Det. w P =4 Det. w P =8 t -5.0 × -4 -2.5 × -4 × -4 × -4 × -4 × -3 C(t)
ParticleStoch. hybridDet. (w P =4) Det. x10
Figure 10: The velocity autocorrelation function for a rigid piston of mas
M/m = 1000 at thermal equilibriumat k B T = 1 (thus, the expected initial value is C (0) = 10 − ), and macro cells of size 4 ×
40 micro cells. Thedeterministic hybrid gives much smaller effective temperature T eff of the piston and a negative dip in C ( t )at short times, but when magnified by an order of magnitude it reveals the correct shape at longer times.The inset focuses on the initial decay of the VACF and shows that even increasing the width of the particleregion to w P = 8 (the whole box has a length of 60 continuum cells) does not help the deterministic hybridmuch whereas the stochastic hybrid gives the correct initial decay. VACF is rather complex due to the interplay of short-time kinetic effects, dissipation, and soundreflections from the walls, however, our focus here is to simply compare against the purely particlesimulations and not to understand all the features in the VACF. The results, shown in Fig. 10,reveal that the piston does not equilibrate at the correct effective temperature T eff = M C (0) /k B in the deterministic hybrid calculations. Notably, just as we found for the massive bead in SectionIV C, the piston has a kinetic energy that is markedly lower (half) than the correct value C (0) = k B T /M = 10 − when fluctuations are not consistently included in the continuum region. Sincethe quasi-equilibrium temperature of the piston plays a crucial role in all of the kinetic theories[68, 69, 71] for the effective heat transfer, it is not surprising that the deterministic hybrid gives thewrong answer. What is even more striking is that increasing the width of the particle subdomain w P rate toward true equilibrium, however, a direct non-equilibrium calculation as we perform in Fig. 9 is more efficientand illustrative for our purpose. It is predicted that the piston equilibrates at a temperature that is approximately the geometric mean of the leftand right temperatures [73]. w P = 2),barely improves the accuracy of the VACF at short times, showing that the whole spectrum offluctuations affects the initial rate of decay of the piston. At the same time, the deterministichybrid does give the correct shape of the VACF at longer times, as revealed by magnifying theVACF for w P = 4 by an arbitrary factor of 10 to bring it in close agreement with the particle resultat longer times . This is also not unexpected since the VACF at longer times is dominated bymechanical vibrations and dissipation, and is analogous to what we find for the dynamic structurefactor when it is calculated by a deterministic hybrid scheme. V. CONCLUSIONS
We have described a hybrid particle-continuum algorithm for simulating complex flows and ap-plied it to several non-trivial problems. The algorithm couples an ideal stochastic particle fluid al-gorithm with a fluctuating hydrodynamic continuum solver using a direct dynamic coupling wherethe continuum solver supplies Dirichlet-like (state) boundary conditions for the particle region,while the particle region supplies von Neumann-like (flux) boundary conditions to the continuumsolver. The continuum solver computes a provisional solution over the whole domain that thengets replaced with the particle data in the particle region and also gets corrected (refluxed) in thecells bordering the particle subdomain. We described the components necessary to extend previousvariants of this rather general coupling methodology [20, 21, 29] to use a recently-developed variantof the DSMC particle method suitable for dense fluids [18, 61], as well as an explicit conservativecompressible fluid solver that accurately accounts for thermal fluctuations in the Navier-Stokesequations [9]. By turning the fluctuating fluxes off in the continuum solver we can trivially trans-form our stochastic hybrid method to a deterministic hybrid method closer to commonly-usedhybrid schemes.In Section IV we applied our stochastic hybrid method to several challenging problems anddemonstrated that it could obtain the correct answer with significantly fewer particles and thussignificantly less computational effort than a purely particle simulation. We used purely particleruns as a “gold” standard against which we judge the accuracy of the hybrid method, consistentwith using the hybrid method for situations in which a purely continuum description is unable tocapture the full physics and a particle method is necessary in some portion of the physical domain. The magnification factor required to bring the long-time VACF for the deterministic hybrid in agreement with thecorrect result decreases with increasing w P , for example, it is ∼
20 for w P = 2 and ∼ w P = 8. N − , where N is the average number of particles in one cell. In Appendix A wedemonstrated that this is not an artifact of the method but rather of an inherent difficulty instochastic methods where instantaneous values are used to estimate means. In section IV B wecalculated the dynamic structure factor using the hybrid method for a quasi two-dimensionalsituation and a large wavevector k that is obliquely incident to the particle-continuum interface,and confirmed that spontaneous sound and entropic fluctuations are transmitted correctly throughthe interface. We also calculated the dynamic structure factor for a finite system bounded by twoadiabatic walls and observed excellent agreement with theoretical calculations given in AppendixB. In Section IV C we studied the diffusive motion of a large spherical neutrally-buoyant beadsuspended in a fluid of I-DSMC particles in three dimensions by placing a mobile particle subdomainaround the suspended bead. We computed the velocity autocorrelation function (VACF) for twobead sizes using the stochastic and deterministic hybrids as well as purely particle simulations.We found that the stochastic hybrid correctly reproduces the VACF over all time scales, while thedeterministic hybrid under-estimates both the kinetic energy of the bead and the magnitude of thetail of the VACF for sufficiently large beads. Finally, in Section IV D we applied the hybrid schemeto a non-equilibrium quasi one-dimensional version of the adiabatic piston problem [67, 68], a classicexample of the importance of fluctuations in establishing global thermal equilibrium. The two-dimensional particle region was placed around the piston and an event-driven algorithm was usedto handle the interaction of the fluid with the piston. We again found that the stochastic hybridwas able to reproduce the purely particle results correctly, while the deterministic calculationsunder-estimated the relaxation substantially for sufficiently massive rigid pistons.Our results for both the VACF and the adiabatic piston clearly demonstrated that a largemassive suspended body cannot equilibrate at the correct Boltzmann distribution unless thermalfluctuations are consistently included in the full domain, including the continuum region in hybridmethods, even if a large particle subdomain is used. This points to an increased importance ofthe long-wavelength, and thus slowly-decaying, hydrodynamic fluctuations. Massive and large3suspended bodies have longer relaxation times and thus it is not surprising that slower-decayingfluctuations play a more prominent role for them than for smaller suspended particles. However,a better theoretical understanding of these observations is necessary in order to establish theimportance of hydrodynamic fluctuations in general. At the same time, our results make it clearthat fluctuations should be included in the continuum region in hybrid methods, consistent withthe particle dynamics, rather than treating the fluctuations as “noise” from which the continuumsolver ought to be shielded.Fluctuating hydrodynamics has successfully accounted for thermal fluctuations in a variety ofproblems. At the same time, however, the nonlinear stochastic partial differential Landau-LifshitzNavier-Stokes equations are mathematically ill-defined and require a cutoff length and/or timescale to be interpreted in a reasonable sense. The linearized equations can be given a well-definedinterpretation, however, they are physically unsatisfactory in that they require a base state to lin-earize around which may itself be time-dependent and unknown or even be affected by fluctuations,as in the adiabatic piston or diffusing shock problems [6, 19]. We have numerically observed thatusing small continuum cells leads to worse results even if the linearized LLNS equations are used,thus formally avoiding the difficulties with the increased relative magnitude of the fluctuations.This is consistent with the expectation that a continuum description is only applicable on lengthscales and time scales sufficiently larger than the molecular size and molecular collision time. Inour experience using more than 75 particles per cell leads to a good match between the hybrid andparticle runs; however, a better theoretical understanding of the proper inclusion of fluctuations inhydrodynamics is a necessary future development.Our implementation is at present serial and our runs are therefore limited by the CPU-intensivecollision procedure in the I-DSMC particle algorithm. Some of the examples we presented utilizedthe mixed event- and time-driven particle algorithm developed in Ref. [30]. The event-drivencomponent of this algorithm is notoriously difficult to parallelize. However, a purely time-drivenparticle algorithm can easily be parallelized, as can the purely continuum solver. We plan toimplement a parallel hybrid scheme in the future in order to enable the study of realistic systemsizes. At the same time, however, reaching long time scales will necessarily require time stepsbeyond the small ones required by particle methods and explicit continuum schemes.4 Appendix A: FINITE-SIZE MISMATCH BETWEEN PARTICLE AND CONTINUUMDESCRIPTIONS
Consider a simple particle/continuum hybrid consisting of a one-dimensional system with anideal particle fluid on one side and a fluctuating continuum solver on the other. At equilibrium themean density and temperature are ρ and T , respectively; the mean fluid velocity is taken as zero.The problem we want to consider is whether the one-sided fluxes of mass, momentum, and energydue to the reservoir particles that enter the continuum subdomain are correct, assuming that themeans and variances of the hydrodynamic variables in the reservoir cells are correct.In the continuum calculation instantaneous density and temperature, ρ and T , are computedand used to generate particles for injection. The one-sided fluxes of mass, momentum, and energyare known from kinetic theory to be F ρ ( ρ, T ) = (cid:114) k πm ρT / F j ( ρ, T ) = k m ρTF E ( ρ, T ) = (cid:114) k m π ρT / . The correct mean equilibrium mass flux is F ρ ( ρ , T ); however, the mean value of the particle fluxdoes not equal this correct value since (cid:104) ρT / (cid:105) = (cid:104) ρ (cid:105)(cid:104) T / (cid:105) (cid:54) = (cid:104) ρ (cid:105)(cid:104) T (cid:105) / ;note that density and temperature are instantaneously uncorrelated. Similarly, the mean energyflux is also incorrect; interestingly the momentum flux is correct since it has a linear dependenceon temperature.To find the errors in the fluxes, we write T = T + δT and for a given power exponent a we have (cid:104) T a (cid:105) = (cid:104) T a (cid:105)(cid:104) (1 + δT /T ) a (cid:105) ≈ (cid:104) T a (cid:105) (cid:20) a ( a − (cid:104) δT (cid:105) T (cid:21) = (cid:104) T a (cid:105) (cid:20) a ( a − γ − N (cid:21) , where N = ρ V c /m is the number of particles in a continuum cell (of volume V c ). From this wehave that, for a monatomic gas ( γ = 5 / (cid:104) F ρ ( ρ, T ) (cid:105) = F ρ ( ρ , T ) (cid:18) − N (cid:19) (cid:104) F J ( ρ, T ) (cid:105) = F J ( ρ , T ) (cid:104) F E ( ρ, T ) (cid:105) = F E ( ρ , T ) (cid:18) N (cid:19) . N − , where N is the number of particlesper macro cell, and therefore the mismatch gets worse as the macroscopic cells become smaller andthe (relative) fluctuations become larger. As our numerical results show, because of this mismatchbetween the particle and hydrodynamic descriptions it will be impossible for the particle andcontinuum regions to reach a common equilibrium state. Appendix B: DYNAMIC STRUCTURE FACTOR WITH ADIABATIC WALLS
Dynamic structure factors are easily calculated in the bulk by using the spatio-temporal Fouriertransform of the linearized LLNS equations [53]. For finite domains, such as slab geometries, thereare results in the literature but they are often restricted to some simplified models or complexnon-equilibrium situations [74, 75]. We therefore derive here the equilibrium dynamic structurefactors for a fluid in-between two adiabatic hard walls by solving the LLNS equations with theappropriate boundary conditions.Boundary conditions change the Hilbert space in which a solution is to be sought and the cor-responding basis functions (eigenfunctions of the generator with the specified BCs). For the LLNSequations with adiabatic boundaries (i.e., slip insulating walls) along the x axes, the appropriatebasis functions are cos( kx ) for density and temperature, where k = pπ/L is the wavevector enu-merated by the wave index p ∈ Z + , and sin( kx ) for the velocities [74, 75], as compared to the onesfor “bulk” conditions (periodic boundaries), exp(2 qπx/L ), with k = 2 qπ/L and wave index q ∈ Z [9, 53]. For thermal walls (constant-temperature stick boundaries) there does not appear to be asimple basis.White noise has a trivial expansion in either the sine or cosine basis sets, namely, all of thecoefficients are i.i.d. Gaussian random variables with mean zero and variance 2. The generator ofthe Navier-Stokes equations separates wavevectors/frequencies into the the same ( k, ω ) equationsas for “bulk” (periodic BCs), and therefore the dynamic structure factor, if expressed in the givenbasis set, has the same familiar form [74, 75]. In particular, ρ ( x, t ) = ρ + ∞ (cid:88) p =1 ρ p ( t ) cos( pπx/L ) = ρ + ∞ (cid:88) p =1 (cid:90) ∞−∞ dωe − iωt ρ p,ω cos( pπx/L ) , where the different p ’s and ω ’s are uncorrelated (cid:10) ρ p,ω ρ (cid:63)p (cid:48) ,ω (cid:48) (cid:11) = 2 ˜ S p,ω δ p,p (cid:48) δ (cid:0) ω − ω (cid:48) (cid:1) , S p,ω = ˜ S k = pπ/L,ω denotes the usual bulk dynamic structure factor.Now, we need to convert this result to the more usual Fourier basis, exp(2 qπx/L ), since this ishow dynamic structure factors are defined. From the Fourier inversion formula, and the orthogo-nality of the cosine basis functions, we have ρ q = 1 L (cid:90) L dxρ ( x, t ) e − qπx/L = 12 ρ p =2 q − iL ∞ (cid:88) p =1 ρ p (cid:90) L dx cos( pπx/L ) sin(2 qπx/L ) . By performing the integration explicitly we get ρ q = 12 ρ p =2 q + 4 iqπ (cid:88) p odd ρ p p − q , giving the dynamic structure factor S ( ρ ) q,ω ≡ S ( ρ ) k =2 πq/L,ω = (cid:10) ρ q,ω ρ (cid:63)q,ω (cid:11) = 12 ˜ S ( ρ ) k,ω + 32 q π (cid:88) p odd ˜ S ( ρ ) p,ω ( p − q ) . (B1)Each of the terms under the sum gives an additional peak at ω = ck = pcπ/L , where p is odd,which arises physically because of the standing waves that appear due to reflections of the soundwaves from the walls. Only the first few terms need to be kept to get most of the power (thetotal integral over ω is one), and it can be shown (by simple numerical comparison or explicitsummation) that the above is equivalent to the more opaque Eq. (12) in Ref. [75] (set γ = 0 forno shear).Equations identical to (B1) hold for temperature and the velocity components parallel to thewall. For the perpendicular velocity ( v x ) a sine basis is more appropriate, and a similar calculationgives the dynamic structure factor for the finite system in terms of the bulk one, S ( v ⊥ ) k =2 πq/L,ω = 12 ˜ S ( v ⊥ ) k,ω + 8 π (cid:88) p odd p ˜ S ( v ⊥ ) p,ω ( p − q ) . Appendix C: HANDLING OF ADIABATIC AND THERMAL WALLS IN THECONTINUUM SOLVER
Solving the LLNS equations with non-periodic boundaries requires some special handling of thestochastic fluxes at the boundaries, which are assumed to coincide with faces of the continuumgrid. As discussed in Refs. [9, 76], the numerical discretization of the Laplacian operator L , thedivergence operator D , and the gradient operator G should satisfy a discrete fluctuation-dissipationbalance condition L = DCG = − DCD ∗ , where C is a dimensionless covariance matrix for thestochastic fluxes that are generated using a random number generator on each face of the grid. For7one-dimension with periodic boundaries it is well known that the standard face-to-cell divergence,cell-to-face gradient, and three-point Laplacian second-order stencils satisfy L = DG and thus C = I (the identity matrix) works and it is in fact trivial to implement algorithmically [9, 76].When boundaries are present, the stencils near the boundaries are modified to take into accountthe boundary conditions.Algorithmically, ghost cells extending beyond the boundaries are used to implement modifiedfinite-difference stencils near the boundaries. The numerical scheme continues to use standarddivergence D (face-to-cell) and gradient (cell-to-face) G = − D (cid:63) stencils but implements a modifiedLaplacian operator due to special handling of the ghost cells. If a Dirichlet condition is imposed ona given variable (e.g., a fixed wall temperature), then the ghost cell value is obtained by a linearextrapolation of the value in the neighboring interior cell (inverse reflection). If a von Neumanncondition is imposed, on the other hand, then the ghost cell value is set equal to the value in theneighboring interior cell (reflection). This gives discrete operators that can be represented by thefollowing banded matrices near the left corner (first cell) of a one-dimensional domain D = − · · · − . . . − · · · ... ... ... . . . , L = − (2 − α ) 1 0 · · · − . . . − · · · ... ... ... . . . , where α = − α = 1 for a von Neumann condition. It is easy to verifythat L = − DCD ∗ is satisfied with the following diagonal scaling matrix C = β · · · . . . · · · ... ... ... . . . , where β = 1 − α. This direct computation shows that in order to satisfy the discrete fluctuation-dissipation bal-ance condition the diagonal element of C corresponding to the cell face that touches the boundaryought to be set to 2 for a Dirichlet and to 0 for a von Neumann condition. This means thatthe corresponding component of the stochastic flux needs to be generated using a random normalvariate of variance 2 for Dirichlet, and set to zero for a von Neumann condition.Finally, for density, the ghost cell values are generated so that the pressure in the ghost cellsis equal to the pressure in the neighboring interior cell, which ensures that there is no unphysical8pressure gradient in the momentum equation across the interface. There is also no stochastic massflux through faces on the boundary independent of the type of boundary condition at the wall. [1] H. Noguchi, N. Kikuchi, and G. Gompper. Particle-based mesoscale hydrodynamic techniques. Euro-physics Letters , 78:10005, 2007.[2] R. Delgado-Buscalioni, K. Kremer, and M. Praprotnik. Concurrent triple-scale simulation of molecularliquids.
J. Chem. Phys. , 128:114110, 2008.[3] G. Hu and D. Li. Multiscale phenomena in microfluidics and nanofluidics.
Chemical EngineeringScience , 62(13):3443–3454, 2007.[4] T. M. Squires and S. R. Quake. Microfluidics: Fluid physics at the nanoliter scale.
Rev. Mod. Phys. ,77(3):977, 2005.[5] N. G. Hadjiconstantinou. The limits of Navier-Stokes theory and kinetic extensions for describingsmall-scale gaseous hydrodynamics.
Physics of Fluids , 18(11):111301, 2006.[6] J. B. Bell, A. Garcia, and S. A. Williams. Numerical Methods for the Stochastic Landau-LifshitzNavier-Stokes Equations.
Phys. Rev. E , 76:016708, 2007.[7] G. De Fabritiis, M. Serrano, R. Delgado-Buscalioni, and P. V. Coveney. Fluctuating hydrodynamicmodeling of fluids at the nanoscale.
Phys. Rev. E , 75(2):026307, 2007.[8] N. K. Voulgarakis and J.-W. Chu. Bridging fluctuating hydrodynamics and molecular dynamics simu-lations of fluids.
J. Chem. Phys. , 130(13):134111, 2009.[9] A. Donev, E. Vanden-Eijnden, A. L. Garcia, and J. B. Bell. On the Accuracy of Explicit Finite-VolumeSchemes for Fluctuating Hydrodynamics. Preprint, arXiv:0906.2425 , 2009.[10] K. Kadau, C. Rosenblatt, J. L. Barber, T. C. Germann, Z. Huang, P. Carles, and B. J. Alder. Theimportance of fluctuations in fluid mixing.
PNAS , 104(19):7741–7745, 2007.[11] O. B. Usta, A. J. C. Ladd, and J. E. Butler. Lattice-Boltzmann simulations of the dynamics of polymersolutions in periodic and confined geometries.
J. Chem. Phys. , 122(9):094902, 2005.[12] P. J. Atzberger, P. R. Kramer, and C. S. Peskin. A stochastic immersed boundary method for fluid-structure dynamics at microscopic length scales.
J. Comp. Phys. , 224:1255–1292, 2007.[13] C. Aust, M. Kroger, and S. Hess. Structure and dynamics of dilute polymer solutions under shear flowvia nonequilibrium molecular dynamics.
Macromolecules , 32(17):5660–5672, 1999.[14] F. J. Alexander and A. L. Garcia. The Direct Simulation Monte Carlo Method.
Computers in Physics ,11(6):588–593, 1997.[15] F. Xijunand N. Phan-Thien, S. Chen, X. Wu, and T. Y. Ng. Simulating flow of DNA suspension usingdissipative particle dynamics.
Physics of Fluids , 18(6):063102, 2006.[16] M. Ripoll, K. Mussawisade, R. G. Winkler, and G. Gompper. Low-Reynolds-number hydrodynamicsof complex fluids by multi-particle-collision dynamics.
Europhys. Lett. , 68(1):106, 2004. [17] S. H. Lee and R. Kapral. Mesoscopic description of solvent effects on polymer dynamics. J. Chem.Phys. , 124(21):214901, 2006.[18] A. Donev, A. L. Garcia, and B. J. Alder. A Thermodynamically-Consistent Non-Ideal StochasticHard-Sphere Fluid. Preprint, arXiv:0908.0510 , 2009.[19] J. B. Bell, J. Foo, and A. L. Garcia. Algorithm refinement for the stochastic Burgers equation.
J.Comput. Phys. , 223(1):451–468, 2007.[20] A. L. Garcia, J. B. Bell, W. Y. Crutchfield, and B. J. Alder. Adaptive Mesh and Algorithm Refinementusing Direct Simulation Monte Carlo.
J. Comp. Phys. , 154:134–155, 1999.[21] S. A. Williams, J. B. Bell, and A. L. Garcia. Algorithm Refinement for Fluctuating Hydrodynamics.
SIAM Multiscale Modeling and Simulation , 6:1256–1280, 2008.[22] M. Moseler and U. Landman. Formation, stability, and breakup of nanojets.
Science , 289(5482):1165–1169, 2000.[23] J. Eggers. Dynamics of liquid nanojets.
Phys. Rev. Lett. , 89(8):084502, 2002.[24] C. Van den Broeck, P. Meurs, and R. Kawai. From maxwell demon to brownian motor.
New Journalof Physics , 7:10, 2005.[25] M. Wu, G. Ahlers, and D.S. Cannell. Thermally induced fluctuations below the onset of Rayleigh-B´enard convection.
Phys. Rev. Lett. , 75(9):1743–1746, 1995.[26] I. Bena, M. Malek Mansour, and F. Baras. Hydrodynamic fluctuations in the Kolmogorov flow: Linearregime.
Phys. Rev. E , 59(5):5503–5510, 1999.[27] I. Bena, F. Baras, and M. Malek Mansour. Hydrodynamic fluctuations in the Kolmogorov flow: Non-linear regime.
Phys. Rev. E , 62(5):6560–6570, 2000.[28] A. Lemarchand and B. Nowakowski. Fluctuation-induced and nonequilibrium-induced bifurcations ina thermochemical system.
Molecular Simulation , 30(11-12):773–780, 2004.[29] S. Wijesinghe, R. Hornung, A. L. Garcia, and N. Hadjiconstantinou. Three-dimensional HybridContinuum-Atomistic Simulations for Multiscale Hydrodynamics.
Journal of Fluids Engineering ,126:768–777, 2004.[30] A. Donev, A. L. Garcia, and B. J. Alder. Stochastic Event-Driven Molecular Dynamics.
J. Comp.Phys. , 227(4):2644–2665, 2008.[31] B. J. Alder and T. E. Wainwright. Studies in molecular dynamics. I. General method.
J. Chem. Phys. ,31:459–466, 1959.[32] B. Duenweg and A. J. C. Ladd. Lattice Boltzmann simulations of soft matter systems.
ArXiv e-prints ,803, 2008.[33] K. T. Mashiyama and H. Mori. Origin of the Landau-Lifshitz hydrodynamic fluctuations in nonequilib-rium systems and a new method for reducing the Boltzmann equation.
J. Stat. Phys. , 18(4):385–407,1978.[34] P. Espa˜nol. Stochastic differential equations for non-linear hydrodynamics.
Physica A , 248(1-2):77–96,1998. [35] L.D. Landau and E.M. Lifshitz. Fluid Mechanics , volume 6 of
Course of Theoretical Physics . Pergamon,1959.[36] M. J. Berger and P. Colella. Local adaptive mesh refinement for shock hydrodynamics.
J. Comp. Phys. ,82:64–84, 1989.[37] J. Bell, M. Berger, J. Saltzman, and M. Welcome. Three-Dimensional Adaptive Mesh Refinement forHyperbolic Conservation Laws.
SIAM Journal on Scientific Computing , 15(1):127–138, 1994.[38] W. Ren. Analytical and numerical study of coupled atomistic-continuum methods for fluids.
J. Comp.Phys. , 227(2):1353–1371, 2007.[39] H. S. Wijesinghe and N. G. Hadjiconstantinou. Discussion of hybrid atomistic-continuum methods formultiscale hydrodynamics.
International Journal for Multiscale Computational Engineering , 2(2):189–202, 2004.[40] W. Ren and W. E. Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics.
J. Comp. Phys. , 204(1):1–26, 2005.[41] T. E. Schwartzentruber and I. D. Boyd. A hybrid particle-continuum method applied to shock waves.
J. Comp. Phys. , 215(2):402–416, 2006.[42] X. B. Nie, S. Y. Chen, W. N. E, and M. O. Robbins. A continuum and molecular dynamics hybridmethod for micro-and nano-fluid flow.
J. Fluid Mech. , 500:55–64, 2004.[43] G. Giupponi, G. De Fabritiis, and P. V. Coveney. Hybrid method coupling fluctuating hydrodynamicsand molecular dynamics for the simulation of macromolecules.
J. Chem. Phys. , 126(15):154903, 2007.[44] R. Delgado-Buscalioni and G. De Fabritiis. Embedding molecular dynamics within fluctuating hydro-dynamics in multiscale simulations of liquids.
Phys. Rev. E , 76(3):036709, 2007.[45] J. Liu, S. Y. Chen, X. B. Nie, and M. O. Robbins. A continuum–atomistic simulation of heat transferin micro-and nano-flows.
J. Comp. Phys. , 227(1):279–291, 2007.[46] S. Chapman and T. G. Cowling.
The Mathematical Theory of Non-Uniform Gases . Cambridge Univ.Press, 1970.[47] A. L. Garcia and B. J. Alder. Generation of the Chapman-Enskog distribution.
J. Comp. Phys. ,140(1):66–70, 1998.[48] S. Chou. Kinetic Flux Vector Splitting for the Navier Stokes Equations.
J. Comp. Phys. , 130:217–230,1997.[49] C. P. Lowe and D. Frenkel. The super long-time decay of velocity fluctuations in a two-dimensionalfluid.
Physica A , 220(3-4), 1995.[50] F. Alexander, A. L. Garcia, and B. J. Alder. Cell Size Dependence of Transport Coefficients in StochasticParticle Algorithms.
Phys. Fluids , 10:1540–1542, 1998. Erratum: Phys. Fluids, 12:731-731 (2000).[51] N. G. Hadjiconstantinou. Analysis of discretization in the direct simulation Monte Carlo.
Physics ofFluids , 12:2634, 2000.[52] A. L. Garcia. Estimating hydrodynamic quantities in the presence of microscopic fluctuations.
Com-munications in Applied Mathematics and Computational Science , 1:53–78, 2006. [53] J. M. O. De Zarate and J. V. Sengers. Hydrodynamic fluctuations in fluids and fluid mixtures . ElsevierScience Ltd, 2006.[54] B. J. Alder and T. E. Wainwright. Decay of the velocity autocorrelation function.
Phys. Rev. A ,1(1):18–21, 1970.[55] M. W. Heemels, M. H. J. Hagen, and C. P. Lowe. Simulating Solid Colloidal Particles Using theLattice-Boltzmann Method.
J. Comp. Phys. , 164:48–61, 2000.[56] T. Ihle and D. M. Kroll. Stochastic rotation dynamics. II. Transport coefficients, numerics, and long-time tails.
Phys. Rev. E , 67(6):066706, 2003.[57] N. Sharma and N. A. Patankar. Direct numerical simulation of the Brownian motion of particles byusing fluctuating hydrodynamic equations.
J. Comput. Phys. , 201:466–486, 2004.[58] J. T. Padding and A. A. Louis. Hydrodynamic interactions and Brownian forces in colloidal suspensions:Coarse-graining over time and length scales.
Phys. Rev. E , 74(3):031402, 2006.[59] P. J. Atzberger. Velocity correlations of a thermally fluctuating Brownian particle: A novel model ofthe hydrodynamic coupling.
Physics Letters A , 351(4-5):225–230, 2006.[60] T. Iwashita, Y. Nakayama, and R. Yamamoto. A Numerical Model for Brownian Particles Fluctuatingin Incompressible Fluids.
Journal of the Physical Society of Japan , 77(7):074007, 2008.[61] A. Donev, A. L. Garcia, and B. J. Alder. Stochastic Hard-Sphere Dynamics for Hydrodynamics ofNon-Ideal Fluids.
Phys. Rev. Lett , 101:075902, 2008.[62] R. Kubo. The fluctuation-dissipation theorem.
Reports on Progress in Physics , 29(1):255–284, 1966.[63] E. H. Hauge and A. Martin-Lof. Fluctuating hydrodynamics and Brownian motion.
J. Stat. Phys. ,7(3):259–281, 1973.[64] R. Zwanzig and M. Bixon. Compressibility effects in the hydrodynamic theory of Brownian motion.
J.Fluid Mech. , 69:21–25, 1975.[65] T. V. Lokotosh and N. P. Malomuzh. Lagrange theory of thermal hydrodynamic fluctuations andcollective diffusion in liquids.
Physica A , 286(3-4):474–488, 2000.[66] T. Iwashita, Y. Nakayama, and R. Yamamoto. Velocity Autocorrelation Function of Fluctuating Par-ticles in Incompressible Fluids: Toward Direct Numerical Simulation of Particle Dispersions.
Progressof Theoretical Physics Supplement , 178:86–91, 2009.[67] E. Lieb. Some problems in statistical mechanics that I would like to see solved.
Physica A , 263(1-4):491–499, 1999.[68] E. Kestemont, C. Van den Broeck, and M. Malek Mansour. The ’adiabatic’ piston: And yet it moves.
Europhysics Letters (EPL) , 49(2):143–149, 2000.[69] J. A. White, F. L. Roman, A. Gonzalez, and S. Velasco. The “adiabatic” piston atequilibrium: Spectral analysis and time-correlation function.
EPL (Europhysics Letters) , 59(4):479–485, 2002.[70] M. Malek Mansour, A. L. Garcia, and F. Baras. Hydrodynamic description of the adiabatic piston.
Phys. Rev. E , 73(1):016121, 2006. [71] M. Cencini, L. Palatella, S. Pigolotti, and A. Vulpiani. Macroscopic equations for the adiabatic piston. Phys. Rev. E , 76(5):051103, 2007.[72] C. Bustamante, J. Liphardt, and F. Ritort. The nonequilibrium thermodynamics of small systems.
Physics Today , 58(7):43–48, 2005.[73] C. Van den Broeck, E. Kestemont, and M. M. Mansour. Heat conductivity by a shared piston.
Euro-physics Letters , 56(6):771–777, 2001.[74] A. L. Garcia, M. Malek Mansour, G. C. Lie, M. Mareschal, and E. Clementi. Hydrodynamic fluctuationsin a dilute gas under shear.
Phys. Rev. A , 36(9):4348–4355, 1987.[75] M. M. Mansour, A. L. Garcia, J. W. Turner, and M. Mareschal. On the scattering function of simplefluids in finite systems.