PaperlayerView exact revision on arXiv

Source-linked AI summary

Staggered Schemes for Fluctuating Hydrodynamics

F. B. Balboa, J. B. BelL, R. Delgado-Buscalioni, A. Donev, T. G. Fai, B. E. Griffith, C. S. Peskin

arXiv:1108.5188v3physics.flu-dyncs.NAmath.NA

TL;DR

Fluctuating hydrodynamics requires numerical schemes that reproduce the covariances of stochastic fields across scales. This paper develops staggered-grid compressible and incompressible schemes with fluctuation-dissipation-preserving spatial discretizations and distinct time integrators. The reported simulations reproduce steady-state covariances and agree with theory for quasi-periodic systems.

  • Problem

    Fluctuating hydrodynamics spans multiple scales and includes fast sound and momentum-diffusion dynamics, making accurate numerical reproduction of fluctuations challenging.

    Source: §II · p. 9

  • Method

    The schemes pair conservative staggered-grid fluxes and discrete fluctuation-dissipation balance with explicit compressible or semi-implicit incompressible time integration.

    Source: §I · pp. 3–4, §VI · p. 43

  • Results

    The schemes reproduce steady-state fluctuation covariances, and quasi-periodic simulation results agree excellently with theory after correcting the transverse wavenumber for spatial discretization artifacts.

    Source: B. Simulations and Results · p. 41, §VI · p. 43

  • Takeaways & Limitations

    The staggered discretization preserves the steady-state fluctuation spectrum while the incompressible method can use large time steps that under-resolve fast grid-scale momentum diffusion.

    Source: §I · p. 3

    For grids much larger than molecular scales, explicit compressible solvers face severe stability restrictions, and zero-Mach variable-density methods remain future work.

    Source: A staggered grid · p. 44

Abstract

from arXiv · show

We develop numerical schemes for solving the isothermal compressible and incompressible equations of fluctuating hydrodynamics on a grid with staggered momenta. We develop a second-order accurate spatial discretization of the diffusive, advective and stochastic fluxes that satisfies a discrete fluctuation-dissipation balance, and construct temporal discretizations that are at least second-order accurate in time deterministically and in a weak sense. Specifically, the methods reproduce the correct equilibrium covariances of the fluctuating fields to third (compressible) and second (incompressible) order in the time step, as we verify numerically. We apply our techniques to model recent experimental measurements of giant fluctuations in diffusively mixing fluids in a micro-gravity environment [A. Vailati et. al., Nature Communications 2:290, 2011]. Numerical results for the static spectrum of non-equilibrium concentration fluctuations are in excellent agreement between the compressible and incompressible simulations, and in good agreement with experimental results for all measured wavenumbers.

I. INTRODUCTION

Thermal fluctuations can produce long-range correlations and rough diffusive interfaces, motivating numerical schemes that preserve their equilibrium statistics. The paper develops staggered-grid compressible and incompressible schemes, with model simplifications and stochastic forcing tailored to the mixing problem.

I. INTRODUCTION

  • The staggered schemes combine conservative fluxes, skew-symmetric advection, and dual discrete divergence and gradient operators to preserve fluctuation-dissipation balance.Compressible flow uses explicit RK3, while incompressible flow uses a semi-implicit method that can take large time steps while retaining steady-state fluctuation covariances.Source: §I · p. 3
  • Imposed gradients can amplify thermal fluctuations and create long-range correlations, producing rough interfaces during diffusive mixing.In microgravity, experiments observed fractal diffusive fronts; the paper focuses on dynamically identical fluid mixtures.Source: §I · pp. 3–4

A. Deterministic Hydrodynamic Equations

  • The mixture model treats the species as nearly identical passive tracers, neglects barodiffusion, and assumes small concentration and constant temperature.A specified constant Soret term is retained to represent the imposed temperature gradient.Source: A. Deterministic Hydrodynamic Equations · p. 5

II. FLUCTUATING HYDRODYNAMICS

  • Stochastic stress and concentration fluxes enter as divergences of Gaussian white-noise fluxes, preserving mass and momentum conservation.The nonlinear white-noise equations require small-scale regularization; the grid or coarse-graining length supplies it.Source: §II · pp. 7–8
  • Compressible equations are stiff because sound waves are fast, while concentration diffusion is slower than momentum diffusion, adding further multiscale dynamics.The paper notes that linearized fluctuating equations are cumbersome and cannot capture possible nonlinear feedback of fluctuations on the mean flow.Source: §II · pp. 7, 9

A. Incompressible Equations

The incompressible formulation projects velocity onto divergence-free fields and uses stochastic forcing whose covariance follows fluctuation-dissipation balance. The paper derives steady-state covariances for equilibrium and a periodic approximation to a concentration gradient.

A. Incompressible Equations

  • The velocity projection is orthogonal onto divergence-free fields, and the stochastic stress forcing is set by fluctuation-dissipation balance.The stress covariance is not uniquely fixed because the projection and divergence operators have non-trivial null spaces.Source: A. Incompressible Equations · pp. 9–10
  • The incompressible projection removes the sound mode, leaving d−1 velocity degrees of freedom per mode.The stochastic stress is taken symmetric because the formulation has no intrinsic-spin equation, although molecular liquids can couple linear and spin momentum.Source: A. Incompressible Equations · p. 10

B. Steady-State Covariances

  • For periodic systems, Fourier transformation decouples modes and yields steady-state covariance matrices, also called static structure factors.These covariances characterize the Gaussian solution of the linearized equations.Source: B. Steady-State Covariances · p. 10
  • At equilibrium, different hydrodynamic variables have no cross-correlations, and their variances are independent of wavevector.The variance of a cell-averaged variable follows by dividing its structure factor by the cell volume.Source: B. Steady-State Covariances · pp. 10–11
  • The non-equilibrium analysis uses a periodic approximation that represents a weak mean concentration gradient through an added velocity-gradient term.The gradient breaks translational symmetry in the actual system, so this approximation linearizes around a uniform background state.Source: B. Steady-State Covariances · p. 11

1. Incompressible Velocity Fluctuations

At equilibrium, incompressible velocity fluctuations have a concentration-gradient-independent covariance proportional to the projection operator. With an imposed concentration gradient, concentration fluctuations are enhanced and develop long-range correlations with velocity fluctuations.

1. Incompressible Velocity Fluctuations

  • At equilibrium, velocity covariance is proportional to the projection operator, and fluctuation amplitude is constant across wavenumbers and reduced by one versus compressible equations.The equilibrium velocity covariance is independent of the concentration gradient.Source: §1 · p. 12
  • Incompressibility removes one kinetic degree of freedom, corresponding to the sound mode, and accounts for the reduced velocity fluctuation amplitude.The removed degree of freedom corresponds to one kBT/2 of kinetic energy.Source: §1 · p. 12

2. Nonequilibrium Fluctuations

  • An imposed concentration gradient enhances concentration fluctuations quadratically and produces concentration–velocity correlations proportional to the gradient, with small-wavenumber divergence indicating long-range correlations.The correlated velocity component is parallel to the concentration gradient.Source: §2 · pp. 13–14

III. SPATIO-TEMPORAL DISCRETIZATION

The temporal schemes target deterministic second-order accuracy and weak accuracy for linearized fluctuating-hydrodynamics equations. Their stochastic updates preserve equilibrium covariances, while using different integrators for compressible and incompressible systems.

A. Temporal Discretization

  • Temporal integrators are higher-order accurate in a weak sense only for the linearized equations, where deterministic accuracy and static-covariance accuracy together determine temporal-correlation accuracy.The schemes are designed to reproduce continuum fluctuation statistics in a weak sense.Source: §III · p. 14, A. Temporal Discretization · pp. 15–16

1. Compressible Equations

  • For compressible equations, a three-stage Runge–Kutta scheme uses two random fluxes per step and achieves third-order weak accuracy for linear additive-noise equations.The scheme is strong-stability-preserving and stable in the inviscid limit without slope-limiting.Source: §1 · p. 16
  • For incompressible equations, a semi-implicit predictor-corrector combines Crank–Nicolson diffusion with explicit treatment of other terms, solving velocity and pressure together without spurious boundary modes.Crank–Nicolson balances numerical dissipation and stochastic forcing identically in the linear setting, preserving the fluctuation spectrum even when fast small-wavelength dynamics are under-resolved.Source: §1 · pp. 18–19
  • The incompressible temporal discretization gives correct steady-state covariances for any time step, although it does not reproduce correct dynamics at large time steps.This holds when advection and diffusion are treated semi-implicitly.Source: §1 · p. 20

B. Spatial Discretization

The spatial formulation uses conservative finite-volume discretization and discrete operators for stochastic advection–diffusion. It is framed as a method-of-lines system in the small-time-step limit.

B. Spatial Discretization

  • The finite-volume discretization represents fields as cell averages and enforces strict local conservation through a conservative divergence.Cell-average changes can be expressed as fluxes through the cell faces.Source: B. Spatial Discretization · p. 20
  • For incompressible flow, the velocity equations can be viewed as constrained stochastic advection–diffusion equations for the separate velocity components.The stochastic stress is split into flux vectors corresponding to the x and y velocity components.Source: B. Spatial Discretization · p. 20
  • In the small-time-step limit, the spatial scheme is a method-of-lines discretization using discrete divergence, gradient, advection, and random-flux operators.The concentration field is represented by finite-volume cell values.Source: B. Spatial Discretization · p. 20

1. Discrete Fluctuation-Dissipation Balance

Discrete fluctuation-dissipation balance requires the spatial scheme to reproduce the correct equilibrium covariance. The staggered grid supports the strict incompressibility constraint needed for velocity statistics, while operator relationships and stochastic flux choices enforce the balance.

1. Discrete Fluctuation-Dissipation Balance

  • A spatial discretization satisfies discrete fluctuation-dissipation balance when it gives the correct equilibrium discrete covariance, making sufficiently small-step discrete statistics consistent with the continuum.At larger time steps, covariance mismatch depends on temporal weak-accuracy order.Source: §1 · p. 21
  • For periodic boundaries, uncorrelated unit-variance random fluxes and negative-adjoint divergence and gradient operators satisfy the balance for dissipative and stochastic terms.The required operator relation is G = −D⋆.Source: §1 · p. 21
  • The discrete advection operator must be skew-adjoint and preserve constants when velocity is discretely divergence-free, so advection does not alter equilibrium fluctuation balance.The stated skew-adjoint condition is S⋆ = −S, with S1 = 0 for constant preservation.Source: §1 · p. 22
  • Incompressible equilibrium covariance is constrained by the discrete projection operator, and the divergence and gradient must satisfy the duality condition G = −D⋆.All discretely incompressible velocity modes should have equal strength at equilibrium for sufficiently small time steps.Source: §1 · p. 22
  • Staggered MAC grids place scalars at cell centers and velocities on faces, offering an alternative to collocated grids where strict discrete divergence-free velocities are difficult to enforce.That collocated-grid difficulty is inconsistent with the incompressible equilibrium covariance condition.Source: §2 · p. 23

3. Stochastic Fluxes

Stochastic fluxes are placed on grid faces and shifted velocity grids, with stress components assigned to grid locations according to their orientation. For concentration-dependent noise, the cell-centered concentration is interpolated to faces by arithmetic averaging, a choice used without theoretical analysis.

3. Stochastic Fluxes

  • The concentration stochastic flux is represented by independent identically distributed numbers on grid faces.Source: §3 · p. 24
  • The concentration-dependent stochastic flux uses concentration interpolated from cell centers to grid faces by arithmetic averaging, a choice made without theoretical analysis.Source: §3 · p. 24
  • Stochastic momentum fluxes occupy shifted velocity grids; in three dimensions, diagonal stress components are at cell centers and off-diagonal components at grid edges.Source: §3 · p. 24
  • For incompressible equations, uncorrelated zero-mean, unit-variance normal variates can generate the stochastic velocity fluxes and yield the correct equilibrium covariances.Alternatively, the stress realization can be made symmetric and traceless; diagonal components are generated by subtracting their average from independent variance-2 normal variates.Source: §3 · p. 24

4. Advection

The advection discretization is designed to be skew-adjoint, supporting discrete kinetic-energy conservation and fluctuation-dissipation balance. Staggered flux construction also preserves divergence-free advection, while the compressible scheme conserves its mass and momentum variables locally.

4. Advection

  • Skew-adjoint advection conserves discrete kinetic energy for inviscid flow and is a key ingredient in fluctuation-dissipation balance.The balance requires the Gibbs-Boltzmann distribution to remain invariant for stochastic velocity dynamics.Source: §4 · p. 25
  • Face-based concentration advection uses a face concentration and velocity; when the original velocity is discretely divergence-free, the shifted advection operator is skew-adjoint.Source: §4 · pp. 25–26
  • The incompressible temporal discretization uses a discretely divergence-free velocity to advect all variables, with the concentration operator discretizing the centered skew-symmetric form of advection.Source: §4 · pp. 25–26
  • With the Soret effect included, concentration is advected using the effective velocity v_adv = v − χ_ST ∇T.Source: §4 · p. 26
  • The compressible scheme uses conservative variables with momentum densities staggered relative to mass densities, and ensures strict local conservation of density, momentum, and partial mass density.Source: §5 · p. 27
  • The compressible equilibrium covariance balance for density–velocity coupling follows from the duality D = −G⋆ between discrete divergence and gradient operators.The text identifies this duality as important for coupling different fluid variables.Source: §5 · p. 28

6. Boundary Conditions

Physical boundaries require modified discrete operators and stochastic-flux treatment that preserve fluctuation-dissipation balance. The boundary rules depend on the condition imposed, while the described boundary discretization retains second-order accurate decay rates for Laplacian modes.

6. Boundary Conditions

  • Boundary modifications to divergence, gradient, and Laplacian operators must preserve fluctuation-dissipation balance; an unsplit Stokes solver leaves that balance unaffected by the discrete incompressibility constraint.Source: §6 · p. 28
  • For cell-centered quantities, the physical boundary is half a cell from the cell center, so ghost cells implement the finite-difference stencils near the boundary.Source: §6 · p. 29
  • For Neumann concentration boundaries, the stochastic boundary flux is zero; for Dirichlet boundaries, the stochastic boundary-flux variance is doubled.Source: §6 · pp. 29–30
  • The boundary Laplacian has continuum-matching normal modes and second-order accurate decay rates, with pointwise second-order accuracy observed even next to the boundary.Quadratic ghost-cell extrapolation can provide formal second-order local accuracy but complicates stochastic-flux handling.Source: §6 · p. 30

IV. IMPLEMENTATION AND NUMERICAL TESTS

The numerical tests examine fluctuation spectra in uniform periodic flows, where discrete fluctuations depend on background-flow magnitude. Accuracy and stability are assessed across wavenumbers and dimensionless advection and diffusion parameters.

IV. IMPLEMENTATION AND NUMERICAL TESTS

  • In uniform periodic flow, the discrete formulation is not Galilean invariant, so fluctuation covariance depends on the background velocity magnitude.Source: §IV · p. 30
  • Explicit advection requires α ≲ 1 for stability, while explicit compressible-flow diffusion requires β, β_c ≤ 1/(2d).The cell Reynolds number r = α/β = V/(νΔx) measures advection relative to dissipation.Source: §IV · pp. 30–31
  • The tests compare discrete Fourier spectra with continuum theory over all discrete wavenumbers, reporting relative variance error for self-correlations and correlation coefficients for cross-correlations.A perfect scheme has normalized concentration spectrum 1 and concentration–velocity cross-correlation 0; the study quantifies deviations as functions of α and β.Source: §IV · p. 31

A. Incompressible Solver

The incompressible solver is assessed through the equilibrium covariance of discretely divergence-free modes and tests of its temporal accuracy. Its observed covariance error decreases quadratically with the time step, while tested spectral deviations remain small.

A. Incompressible Solver

  • The incompressible scheme is implemented in IBAMR, with SAMRAI managing parallel Cartesian grids and PETSc providing iterative Krylov solvers.The Stokes linear solver accounts for most computational effort; hypre geometric multigrid is used to solve the linear systems efficiently.Source: A. Incompressible Solver · pp. 31–32
  • Equilibrium accuracy is assessed using covariance of an orthonormal basis of discretely divergence-free modes, whose amplitudes should be unity across wavenumbers, including with physical boundaries.For equilibrium, the random coefficients in this basis should be uncorrelated and have unit variance.Source: A. Incompressible Solver · p. 32
  • With r = 1 and α = 0.5 or 0.25, deviation from unity is less than 5% even at the larger time step; covariance is correct at all wavenumbers as α→0.The plotted quantity is the spectral power of the first solenoidal mode as a function of wavenumber.Source: A. Incompressible Solver · p. 32, Figure 2 · p. 33
  • The observed velocity-variance error decreases quadratically with α, and concentration gives identical results, confirming discrete fluctuation-dissipation balance and weak second-order temporal accuracy.The numerical results confirm the predicted quadratic convergence.Source: A. Incompressible Solver · pp. 32–33

B. Compressible Solver

The compressible solver reproduces equilibrium covariances with weak third-order temporal accuracy, and its static and small-wavenumber dynamic structure factors agree with theoretical expectations. The paper also applies staggered schemes to giant fluctuations in microgravity, using approximations to make simulations and theory tractable.

B. Compressible Solver

  • The explicit compressible GPU implementation runs about four times faster than a CPU-based code on 32 AMD cores using MPI.The CUDA implementation assigns one thread per cell and accesses that cell and its neighbors.Source: B. Compressible Solver · p. 33
  • Self-correlations remain close to unity and cross-correlations nearly vanish; density fluctuations have up to 5% relative error at the largest wavenumbers.The equilibrium spectra are examined for a uniform periodic system with imposed background flow.Source: B. Compressible Solver · p. 34, Figure 4 · p. 35
  • Cubic covariance-error convergence verifies discrete fluctuation-dissipation balance and weak third-order temporal accuracy for the compressible scheme.Figure 3 compares covariance errors with the predicted O(∆t^3) convergence.Source: B. Compressible Solver · pp. 34–35

1. Dynamic Correlations

  • At the smallest resolved wavenumber, numerical and theoretical dynamic structure factors agree very well, while larger wavenumbers show discretization artifacts.Peak positions and widths are governed by the effective rather than true wavevector.Source: §1 · pp. 36–37, Figure 5 · p. 36

V. GIANT FLUCTUATIONS

  • The authors perform the first incompressible simulations of diffusive mixing in microgravity, comparing the experimental setting with fluctuating hydrodynamics.The experimental measurements were made aboard a satellite.Source: §V · p. 37
  • To make steady-state spectra feasible to compute, the simulations reduce the Schmidt number while keeping sχ^2 fixed and s ≫ 1.The experimental Schmidt number is approximately 3 · 10^3, making direct temporal simulation infeasible for long averaging runs.Source: §V · p. 38

A. Approximate Theory

  • The quasi-periodic model uses periodic boundaries and an imposed concentration-gradient source, but is used only for testing and theoretical analysis, not experimental comparison.Its mean concentration profile is uniform, unlike the exponential profile in experiments and simulations with physical boundaries.Source: A. Approximate Theory · p. 39
  • No-slip boundaries suppress fluctuations at scales larger than their separation; after the gradient is established, giant fluctuations evolve on a typical timescale of about 1000 s.The concentration snapshots illustrate evolution from an initially fluctuation-free state toward a steady state whose typical length scale is set by the finite domain.Source: A. Approximate Theory · p. 40, Figure 6 · p. 40

B. Simulations and Results

The simulations compare concentration spectra from compressible and incompressible schemes with theory and experimental observations. The conclusions emphasize discrete fluctuation-dissipation balance, the advantages of staggered grids for incompressible flow, and the scope of future modeling.

B. Simulations and Results

  • Explicit concentration boundary conditions with semi-implicit diffusion do not conserve total concentration, but the observed average-concentration change is small for the parameters used.The authors note that implicit boundary handling or explicit concentration diffusion could correct this behavior.Source: B. Simulations and Results · p. 41
  • Both panels show excellent agreement between theory and numerical results for quasi-periodic systems, with ˜k⊥ correcting most spatial-discretization artifacts.Figure 7 compares the steady-state concentration spectrum averaged along the y-axes in two and three dimensions, using the ratio of observed to predicted spectra.Source: B. Simulations and Results · p. 41, Figure 7 · p. 42
  • The semi-implicit incompressible scheme gives the correct static concentration covariance at all time-step sizes, and simulations generally use viscous CFL numbers β = 1 or 2.Some large systems use β = 5; concentration boundary conditions are imposed explicitly, causing a small average-concentration change for the parameters used.Source: B. Simulations and Results · pp. 41–42
  • For non-periodic systems, single-mode Galerkin theory over-predicts concentration fluctuations at smaller wavenumbers in both dimensions, while results for Schmidt numbers 10 and 20 differ little.Compressible and incompressible two-dimensional simulations agree excellently; the figure compares both formulations, boundary conditions, Schmidt numbers, and experimental data.Source: B. Simulations and Results · p. 42
  • Within arbitrary normalization, numerical results agree well with experimental observations across accessible wavenumbers, although the smallest-wavenumber comparison is limited by boundary mismatch and measurement uncertainty.Agreement at small wavenumbers improves as the simulation Schmidt number increases.Source: B. Simulations and Results · pp. 42–43

VI. CONCLUSIONS

  • Discrete fluctuation-dissipation balance uses L = −DD⋆ with boundary corrections, while centered skew-adjoint advection adds no fluctuation forcing or dissipation.The spatial discretization is designed so each mode is equally forced and dissipated, consistent with the continuum balance for small time steps.Source: §VI · p. 43

A staggered grid

  • For incompressible flow, a staggered grid enables a robust idempotent projection that maintains discrete fluctuation-dissipation balance and avoids artificial pressure modes at physical boundaries.On periodic systems, the temporal method is equivalent to a classical Crank-Nicolson-based projection method.Source: A staggered grid · p. 44
  • Zero-Mach methods for variable-density fluctuating hydrodynamics remain a future challenge for modeling buoyancy in Earth-based giant-fluctuation experiments.The compressible solver uses GPUs, and the authors defer analysis of temporal correlations to future work.Source: A staggered grid · p. 44

Appendix A: Implicit Midpoint Rule as a Gibbs Sampler

The appendix frames the implicit midpoint rule as a sampler for the equilibrium Gibbs distribution and analyzes timestep choices for accurate sampling. It also shows that the incompressible MAC projection and Crank-Nicolson discretization reproduce the correct equilibrium velocity covariance at any timestep.

Appendix A: Implicit Midpoint Rule as a Gibbs Sampler

  • For steady-state observables, a sampler need only reproduce the equilibrium Gibbs distribution, not the correct dynamics.The appendix uses this distinction to motivate treating the implicit midpoint or Crank-Nicolson discretization as an MCMC sampler.Source: Appendix A · pp. 45–46
  • The implicit midpoint rule samples the correct equilibrium distribution exactly, and its Metropolis-Hastings acceptance probability is α = 1 for any timestep.The rule can therefore be viewed as an MCMC algorithm whose equilibrium distribution is the target Gibbs distribution.Source: Appendix A · p. 46
  • The timestep can be chosen without biasing equilibrium sampling, but statistical accuracy depends on the observable’s autocorrelation time.For the simulations discussed, β ≈ 2 balances coverage across wavenumbers because the fluctuation spectrum has power-law tails ∼k^-4.Source: Appendix A · pp. 46–47
  • For small dimensionless wavenumbers, increasing the timestep improves sampling, whereas at large wavenumbers it reduces statistical accuracy.The appendix relates this trade-off to Crank-Nicolson being A-stable but not L-stable.Source: Appendix A · p. 47

Appendix B: Fluctuation-Dissipation Balance for Incompressible Flow

  • For incompressible flow, the MAC projection is self-adjoint and idempotent, giving the correct velocity covariance with Crank-Nicolson at any timestep.The covariance is Sv = P.Source: Appendix B · p. 48
  • The incompressible velocity covariance result follows from fluctuation-dissipation balance when the projection P is self-adjoint and idempotent.In the MAC discretization, P = I − D⋆(DD⋆)^−1D meets both conditions.Source: Appendix B · p. 48
Loading 1108.5188v3…