# Projection Method for the Fluctuating Hydrodynamics Equations

Computational fluctuating hydrodynamics aims at understanding the impact of thermal fluctuations on fluid motions at small scales through numerical exploration. These fluctuations are modeled as stochastic flux terms and incorporated into the classical Navier-Stokes equations, which need to be solved numerically. In this paper, we present a novel projection-based method for solving the incompressible fluctuating hydrodynamics equations. By analyzing the equilibrium structure factor spectrum of the velocity field, we investigate how the inherent splitting errors affect the numerical solution of the stochastic partial differential equations in the presence of non-periodic boundary conditions, and how iterative corrections can reduce these errors. Our computational examples demonstrate both the capability of our approach to reproduce correctly stochastic properties of fluids at small scales as well as its potential use in the simulations of multi-physics problems.

## Authors

• 1 publication
• 1 publication
• 1 publication
03/29/2020

### An iterative splitting method for pricing European options under the Heston model

In this paper, we propose an iterative splitting method to solve the par...
06/14/2020

### TURB-Rot. A large database of 3d and 2d snapshots from turbulent rotating flows

We present TURB-Rot, a new open database of 3d and 2d snapshots of turbu...
07/30/2019

### Iterative and Non-iterative Splitting approach of a stochastic Burgers' equation

In this paper we present iterative and noniterative splitting methods, w...
08/15/2021

### Convergence study of IB methods for Stokes equations with non-periodic boundary conditions

Peskin's Immersed Boundary (IB) model and method are among the most popu...
03/24/2021

### Fluid Flow along the Riga Plate with the Influence of Magnetic Force in a Rotating System

The fluid flow along the Riga plate with the influence of magnetic force...
03/21/2020

### Parameter robust preconditioning by congruence for multiple-network poroelasticity

The mechanical behaviour of a poroelastic medium permeated by multiple i...
10/25/2019

### Bayesian identification of a projection-based Reduced Order Model for Computational Fluid Dynamics

In this paper we propose a Bayesian method as a numerical way to correct...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Owing to the perpetual miniaturization of engineered systems, the smallest scale of interest for engineers and scientists has been continuously shrinking, leading to the emergence of nanometric devices and a growing curiosity into the fundamental properties of nanofluidic systems nanofluidicOnTheRise20201 . This trend has brought various exciting perspectives and novel ideas to the scientific community, including bioengineered kinesin-microtubule systems lopes_membrane_2019 , nanofabricated devices for biomolecule applications B917759K , nano heat transport technologies nanoheatreview ; Wang2009 , brain-on-a-chip platforms 10.3389/fbioe.2019.00100 , and micro-sensing devices microsensing , to name but a few. These systems display exotic behaviors and acquire unique features because their characteristic length scales are small, for example, comparable to the Debye length, the size of biomolecules, or even the slip length Abgrall2008 .

In particular, since the characteristic length and time scales of the systems investigated are no longer widely separated from those of the underlying molecular systems, the validity of any deterministic continuum models is questionable and should be investigated. Surprisingly, experimental doi:10.1063/1.449693 ; doi:10.1063/1.465059 ; doi:10.1126/science.1074481 ; PhysRevB.75.115415 ; PhysRevLett.96.086105 ; PhysRevLett.91.166104 , theoretical B909366B , and numerical molhydro ; PhysRevLett.102.184502 ; B616490K ; PhysRevLett.94.026101 studies have found that the Navier–Stokes equations remain largely valid down to a few nanometers. Specifically, the classical hydrodynamics equations can predict the average nanofluidics behavior near the equilibrium or in simple steady states even at the nanometer scale. However, thermal fluctuations appearing in the instantaneous continuum fields are no longer negligible as the system size becomes smaller. When these fluctuations interact with nonlinearity in the system, the entire dynamics of the fluid system may not be correctly captured by the deterministic continuum description. For example, in the giant fluctuation experiment in space GFexp ; GFexp2 , concentration fluctuations have been observed to grow to the macroscale (i.e. sizes ranging up to millimeters and relaxation times as large as hundreds of seconds) under the presence of the concentration gradient due to the coupling with random advection. Therefore, for a complete and accurate representation of nanofluidic and sub-nanofluidic systems, thermal fluctuations must be accounted for.

The use of stochastic partial differential equations (SPDEs) to describe fluid dynamics dates back to Landau and Lifshitz LandauLifshtz1959 . They proposed to incorporate stochastic fluxes to each dissipative process (e.g. momentum diffusion) in the Navier–Stokes equations to correctly model the effects of thermal fluctuations. This fluctuating hydrodynamics (FHD) approach has been successfully applied to describe various phenomena induced by hydrodynamic fluctuations. Until significant progress in the computational FHD approach has been made for the past two decades (see below), most accomplishments of the FHD theory were made by analytical methods OrtizDeZarateSengers2006 . While analytical approaches have provided insightful explanations, they are mostly limited to simple nonequilibrium situations and moreover, they rely on several assumptions (e.g. linearization and periodic boundary conditions). We note that stochastic terms can also be added to the deterministic fluid equations in the context of uncertainty quantification MathelinHussainiZang2005 ; Wan2007 to represent a large variety of noisy contributions (e.g. uncertainty on the boundary conditions PhysRevLett.92.154501 ).

As mentioned above, significant progress in the computational FHD approach has been made for the past two decades. Here we focus on the PDE-based (as opposed to particle-based) approaches for homogeneous fluids. The Landau–Lifshitz Navier–Stokes equations (i.e. compressible FHD equations) were numerically solved for the one-species CCSE_PRE2007 and binary mixture CCSE_ESAIM2010 cases as well as the multi-species case CCSE_PRE2014 , which was extended to include stochastic reactions CCSE_JCP2015 . The incompressible FHD formulation BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 was extended to quasi-incompressible fluids by the low Mach number formulation for the binary mixture CCSE_CAMCOS2014 ; CCSE_CAMCOS2015 and muiti-species cases CCSE_PhysFluids2015 . The quasi-incompressible case was extended to include stochastic reactions CCSE_JChemPhys2018 as well as charges (i.e. electrolyte ions) CCSE_PRF2016 ; CCSE_PRF2019 . To construct and analyze numerical schemes for FHD equations, various advanced deterministic PDE and computational fluid dynamics (CFD) techniques have been applied and extended. Spatial discretization based on the finite-volume approach and the stochastic version of the method of lines were introduced, and the structure factor analysis technique was developed for the systematic construction of stochastic numerical methods DonevVandenEijndenGarciaBell2010 . For (quasi-)incompressible FHD equations, staggered spatial discretization schemes (i.e. using a grid with staggered momenta) were developed BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 ; VoulgarakisChu2009 . In addition, to solve the stochastic Stokes problem, which is a saddle-point linear system, using the generalized minimal residual method (GMRES), an efficient variable-coefficient finite-volume Stokes solver was developed CaiNonakaBellGriffithDonev2014 . Several time integrators for FHD equations (e.g. semi-implicit schemes CCSE_CAMCOS2015 ) were also constructed and analyzed DelongGriffithVandenEijndenDonev2013 . While it is not possible to give a complete summary of applications and extensions of the FHD approach here, we note that the FHD approach has been applied to reaction-diffusion systems Atzberger2010 ; CCSE_JChemPhys2017 and coupled to kinetic Monte Carlo SelmiMitchellManoranjanVoulgarakis2017 and molecular dynamics VoulgarakisChu2009 . Hydrodynamic couplings between microstructures or ions were also considered using the stochastic Eulerian-Lagrangian method Atzberger2011 ; wang_fluctuating_2018 , the boundary integral formulation BaoRachhKeavenyGreengardDonev2018 , and the immersed boundary approach CCSE_PRF2021 .

This paper presents a projection-based method for solving the incompressible FHD equations, with the dual intent to provide an alternative simulation strategy and illustrate how traditional CFD techniques can be adapted to incorporate thermal fluctuations. The projection method CHORIN196712 uses the Hodge decomposition to decouple the fluid equations and update the solution in two steps. First, the momentum contributions are used to advance the velocity field, which is then projected on the divergence-free space to enforce the incompressibility condition and recover the pressure field. Owing to this decoupling, the projection method alleviates the need for constructing and solving a monolithic system containing the coupled hydrodynamic equations by forming and solving two smaller systems. Using traditional data structures and discretization strategies GUITTET2015 , one can guarantee that these systems are symmetric positive definite and therefore can efficiently solve them using classical iterative methods, which can be accelerated using parallel architectures EGAN2021110084 . However, when physical boundary conditions such as the no-slip boundary condition are used, it is well known that decoupling causes errors near the boundaries ELiu1995 ; ELiu2003 . To avoid this issue, (semi-)implicit schemes for solving the incompressible FHD equations have been developed mostly by solving the coupled system CCSE_CAMCOS2015 (note, however, that the projection approach is still employed as a preconditioner CaiNonakaBellGriffithDonev2014 ). Nevertheless, since the computational advantage of using the projection method in FHD simulations is expected to be significant, particularly for large-scale simulations, we propose a projection-based method with iterative boundary corrections and perform a systematic numerical analysis based on the equilibrium structure factor.

The rest of the present paper is organized as follows. We start in section 2 by introducing the stochastic incompressible Stokes equation and its spatial discretization on uniform staggered grids. In section 3, we introduce the steady-state covariance and static structure factor that will be used to analyze our numerical schemes. Our projection method is presented in section 4, analyzed in section 5, and numerically validated in section 6. In section 7, we employ our method to simulate giant fluctuations. We conclude in section 8.

## 2 Stochastic Incompressible Stokes Equation

We introduce here an SPDE for the velocity field of an incompressible fluid and discuss its spatial and temporal discretizations. We consider the equilibrium case, where the SPDE can be linearized. The resulting stochastic incompressible Stokes equation is presented in section 2.1. Its spatial discretization based on the finite-volume approach is described in section 2.2. As an example of a numerical scheme that does not use the projection method approach, a temporal integration scheme based on the Crank–Nicolson approximation is given in section 2.3.

### 2.1 Continuum Equation

The fluctuating behavior of an incompressible fluid in equilibrium can be modeled by the following FHD equations:

 ρ∂u∂t+ρu⋅∇u=μΔu−∇π+√2kBTμ∇⋅Σ, (1a) ∇⋅u=0, (1b)

where and denote respectively the velocity and pressure fields, and are respectively the mass density and temperature of the fluid and are taken as constant, is the Boltzmann constant, and

is the dynamic viscosity. The tensor field

denotes the stochastic momentum flux.

When velocity fluctuations can be assumed to be small, the self-advected term , which is of second order in , can be neglected. Under this assumption, using the kinematic viscosity and the orthogonal projection onto the space of divergence-free velocity fields, we linearize equation (1) as

 ∂u∂t=P[νΔu+√2kBTνρ∇⋅W]. (2)

Here, we assume that

is a spatiotemporal Gaussian white noise tensor field with independent components having covariance

 ⟨Wij(r,t)Wi′j′(r′,t′)⟩=δii′δjj′δ(r−r′)δ(t−t′). (3)

In principle, the symmetrized form (i.e.  should have been used. However, what matters in the Fokker–Planck description is the covariance of the projected stochastic forcing (see (4) below) and the use of the symetrized form is not necessary for incompressible flow with constant viscosity BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 ; DelongGriffithVandenEijndenDonev2013 . We also assume that the initial velocity is divergence-free, i.e. . This assumption implies that remains to be divergence-free for all .

In this paper, we consider the periodic boundary condition and the no-slip boundary condition (i.e.  on the boundary). For both boundary conditions, the divergence operator and gradient operator satisfy

, where star denotes an adjoint of a matrix or linear operator. Hence, the vector Laplacian operator

is self-adjoint, i.e. , and the projection operator with being the identity operator is indeed an orthogonal projection, i.e.  and . In addition, using and , the covariance of the projected stochastic forcing is expressed as BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012

 ⟨(PDW(t))(PDW(t′))∗⟩=−δ(t−t′)PLP. (4)

Therefore, it can be seen in (2) that the stochastic term is linked to the viscous term LandauLifshtz1959 by the fluctuation-dissipation balance FluctuationDissipation .

### 2.2 Spatial Discretization

To spatially discretize the stochastic incompressible Stokes equation (2), we employ the staggered-grid discretization on uniform mesh developed in BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 . While we consider the two-dimensional case in this paper for clarity, the three-dimensional case is essentially the same. We denote the spatially discretized equation as

 dUdt=P[νLU+√2kBTνρΔVDWW], (5)

where and are respectively the discretizations of and the stochastic tensor field ; and are respectively the discrete projection and vector Laplacian operators; is the discrete divergence operator acting on ; and is the volume of a cell. The data layout and discrete operators are illustrated in Figure 1 and described next. In section 2.2.1, we first describe variables and operators, assuming periodic boundary conditions. In section 2.2.2, we explain modifications needed to impose the no-slip boundary condition.

#### 2.2.1 Variables and Operators

According to the marker-and-cell layout HarlowWelch1965 , the velocity components ( and ) are stored at the faces, whereas the pressure and Hodge variable (denoted as ) components are stored at the cell centers. The discrete divergence operator and discrete gradient operator are defined as follows. At a cell center , the divergence of the velocity, , is constructed using the second-order centered finite difference. At face , the Hodge gradient is also computed using the second-order centered finite difference. With these definitions, desired relations that hold in the continuous case are still valid. First, the discrete gradient and divergence operators obey the duality relation . From the definition of the discrete projection operator , it can be easily seen that is indeed an orthogonal projection, i.e.  and . The discrete Laplacian of the velocity at face , denoted by , is computed using the standard second-order 5-point stencil.

The diagonal components of the stochastic stress tensor ( and ) are stored at the cell centers, while the off-diagonal components ( and ) are stored at the nodes of the mesh (see Figure 1). These stochastic terms are constructed as follows. Since each component of the stochastic noise is a distribution, it cannot simply be evaluated at any given point and must be interpreted in the integral form. That is, its discretization is constructed as the spatial average over the volume of size , centered where is stored:

 ¯¯¯¯¯¯Wij(t)=1ΔV∫VijWij(r,t)dr. (6)

Hence, each component

defined at each cell center or node is an independent Gaussian white noise process with variance

. To explicitly express the dependence of the magnitude of fluctuations on , we introduce normalized stochastic processes . The covariance of is expressed as

 ⟨W(t)W(t′)⟩=CWδ(t−t′). (7)

While

is simply the identity matrix in the periodic boundary case, we will see that

needs to be modified for the no-slip boundary condition.

One of the crucial issues for spatial discretization is that the discrete system should reproduce a correct equilibrium distribution that is expected from the continuous case. Since the fluctuation-dissipation balance principle dictates the equilibrium in the continuous case, it is required that its discrete version should be satisfied when a spatial discretization is constructed. In other words, the discrete fluctuation-dissipation balance dictates the choice of . Since is a discretization of a divergence operator, it is natural to base its construct on , keeping in mind that acting on is defined at the cell centers, while acting on is defined at the faces. Using the second-order centered finite differences, the discrete stochastic divergence is constructed. Then it can easily be seen that the following discrete fluctuation-dissipation balance

 L=L∗=−DWCWD∗W (8)

is satisfied. The properties of the discretized operators and are important for computing the steady-state covariance in sections 3 and 5.

#### 2.2.2 Boundary Conditions

In the presence of non-periodic boundaries, the discrete operators defined above need to be modified near the boundaries. For the no-slip boundary condition, the velocity component normal to the boundary is zero at the boundary. Hence, the velocity components at faces on the boundary are set to zero and not included as independent degrees of freedom (see Figure

2). Then, no values in cells outside the physical domain (i.e. ghost cells) are needed to define . For the projection step, this zero normal velocity condition implies that the homogeneous Neumann boundary condition should be chosen for the Hodge variable ELiu1995 . Hence, ghost cell values for the Hodge variable are set to be equal to the values in the neighboring interior cells, and the resulting operator satisfies the duality relation with  BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 . Therefore, and continue to hold for .

To define the discrete vector Laplacian , ghost cells are needed for parallel velocity components that are half a grid spacing away from the boundary (see Figure 2). When the tangential velocity is prescribed at the nodes on the boundary, the corresponding ghost cell value is determined by the linear extrapolation BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 . Then, can be constructed using the 5-point stencil, and can be expressed as

 LU=L0U+B¯U, (9)

where denotes the prescribed tangential velocity. Since we consider the homogeneous Dirichlet boundary condition (i.e. ), becomes , and thus does not appear. However, it is noted that the intermediate velocity appearing in a projection method may not satisfy the prescribed boundary condition. To develop and analyze our numerical schemes, we will use this representation for nonzero tangential velocities. We note that the size of depends on how is represented. For the data layout of , we use the same one as , where each tangential velocity prescribed on the boundaries (empty squares in Figure 2) is stored at the location of the face (the closest filled triangle) half a grid spacing inward from the actual location of the prescribed velocity on the boundaries. In this setting, the size of is identical to that of .

Finally, since is modified, also needs to be modified to satisfy the discrete fluctuation-dissipation balance (8). This can be achieved by setting the variance of stochastic fluxes (affecting the tangential velocity components) at nodes on the boundary to twice that used for the interior fluxes BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 . In other words, is still a diagonal matrix with most diagonal elements being one but diagonal elements corresponding to those nodes on the boundary become two.

### 2.3 Temporal Discretization

Following DelongGriffithVandenEijndenDonev2013 , we discretize (5) in time as

 Un+1=P[Un+νΔtL(Un+Un+12)+√2kBTνΔtρΔVDWWn]. (10)

Here, superscripts denote timesteps. Since the covariance of is proportional to , the collection of the Gaussian white noise processes has been replaced by , where each component of defined at each spatial point at each time step is an independent standard normal random variate. As explained in section 2.2.2, the variance of needs to be doubled on the no-slip boundary to satisfy the discrete fluctuation-dissipation balance (8).

This scheme is obtained from the Crank–Nicolson approximation with the assumption . Note that the resulting also stays in the projected space. As explained in section 3.2, this scheme does not introduce any time discretization error in the equilibrium covariance (or equivalently, the equilibrium structure factor) DelongGriffithVandenEijndenDonev2013 . In order to implement this scheme, however, a linear solver CaiNonakaBellGriffithDonev2014 for the following coupled system is required DelongGriffithVandenEijndenDonev2013 :

 (I−12νΔtL)Un+1+ΔtGΠn+12=(I+12νΔtL)Un+√2kBTνΔtρΔVDWWn, (11a) DUn+1=0. (11b)

The main goal of this paper is to solve (10) not using a linear solver for the saddle-point system (11) but using a projection method.

## 3 Steady-State Covariance and Static Structure Factor

In this section, we introduce quantities that characterize the behavior of fluctuations in the equilibrium Stokes system and corresponding discretized systems. Since the mean velocity is zero at equilibrium, we focus on the second moment. The steady-state covariance measures the covariance of velocities at two physical locations for a system in equilibrium. Since the velocity fields are represented by

, , and for the continuum, spatially discretized, and fully discretized cases, respectively, the corresponding steady-state covariances are defined as

 Continuum R =limt→∞⟨u(t)u∗(t)⟩, (12a) Spatially discretized RΔx =limt→∞⟨U(t)U∗(t)⟩, (12b) Fully discretized RΔx,Δt =limn→∞⟨Un(Un)∗⟩. (12c)

For each case, the static structure factor is defined as the Fourier transform of the steady-state covariance. We note that, for the no-slip boundary case, the velocity field is mirrored and the resulting field in the extended domain is used (see Figure

11) since the velocity field in the original domain is not periodic.

These quantities can be used to investigate the accuracy of spatial discretization and numerical schemes. In section 5, we analytically investigate our projection-method-based schemes by computing the steady-state covariance. In section 6, we numerically investigate those schemes by computing the static structure factor. In this section, we derive the main analytic results for the continuum and spatially discretized cases (i.e.  and ) as well as the Crank–Nicolson scheme (10) (i.e. ). While those results are known BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 , we present them along with derivations, as both aspects are essential for the analysis of our projection-method-based schemes. In A, we develop a systematic procedure to construct and solve linear systems from which the steady-state covariance can be uniquely determined.

### 3.1 Continuum and Spatially Discretized Cases

Using the result of A, the steady-state covariance of the continuum equation (2) is given as the unique solution of

 νPLR+νRLP=−2kBTνρPDD∗P, (13a) PR=RP=R. (13b)

The physical intuition that the equilibrium covariance must be proportional to the projection operator, suggests that

 R=kBTρP. (14)

Using the properties of and mentioned in section 2.1 (i.e. , , , ), it can be easily shown that (14) indeed satisfies (13).

The static structure factor is obtained from the Fourier transform of the steady-state covariance (with normalization DonevVandenEijndenGarciaBell2010 with respect to the volume of the system ):

 S=V⟨ˆu(t)ˆu∗(t)⟩=kBTρˆP, (15)

where a hat denotes the Fourier transform. For the periodic boundary case, the structure factor has a compact form

 S(k)=kBTρ[I−kk∗k⋅k]. (16)

It is easy to see from (16) that all divergence-free modes have the same spectral power at equilibrium in the periodic boundary case. In fact, we note that (15) implies the same conclusion for a general case where aforementioned properties of and hold BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 .

For the spatially discretized equation (5), a similar argument can be made since the discrete operators and are constructed so that they preserve all aforementioned properties of and . Because in this case the fluctuation-dissipation balance is given in the form of (8), the linear system that uniquely determines is

 νPLRΔx+νRΔxLP=−2kBTνρΔVPDWCWD∗WP (17a) PRΔx=RΔxP=RΔx (17b)

and, similarly to the continuum steady-state expression (14), the discrete steady-state covariance is given as

 RΔx=kBTρΔVP. (18)

### 3.2 Crank–Nicolson Scheme

Introducing , the Crank–Nicolson scheme (10) can be written as

 Un+1=A−1−A+Un+√2kBTνΔtρΔVA−1−PDWWn. (19)

Hence, using the result of A, the steady-state covariance is given as the unique solution of

 (A−1−A+)RCNΔx,Δt(A−1−A+)∗−RCNΔx,Δt=−2kBTνΔtρΔV(A−1−PDW)CW(A−1−PDW)∗, (20a) PRCNΔx,Δt=RCNΔx,ΔtP=RCNΔx,Δt. (20b)

By observing that

 A+PA∗+−A−PA∗−=2νΔtPLP, (21)

and using the discrete fluctuation-dissipation balance (8), one can show that

 RCNΔx,Δt=kBTρΔVP (22)

satisfies (20). The results (18) and (22) show that the Crank–Nicolson scheme does not introduce any temporal integration errors to the steady-state covariance.

## 4 Construction of Projection Methods

We present here how our projection methods are constructed to compute the numerical solution of the (spatially discretized) stochastic incompressible Stokes equation (5). The resulting schemes for the periodic boundary and no-slip boundary cases are given in Schemes 1 and 2, respectively. Our projection-method-based schemes solve the Crank–Nicolson update (10) using the operator splitting approach. In section 4.1, we consider the periodic boundary case, where the simple splitting exactly solves (10). In section 4.2, we first discuss issues that arise when the simple splitting is applied to the no-slip case and introduce the idea of iterative boundary corrections. In section 4.3, we present our iterative scheme for the no-slip case.

### 4.1 Periodic Boundary Case

For periodic boundary conditions, the Crank–Nicolson update (10) can be exactly implemented by the simple projection-based time-splitting given in Scheme 1. Since both and obey periodic boundary conditions, also satisfies periodic boundary conditions.

Using the fact that and commute in the periodic boundary case, we can show that Scheme 1 does not introduce any splitting error to solve (10). We first observe that Scheme 1 can be written as

 Un+1=P~U=P[Un+νΔtL(Un+~U2)+√2kBTνΔtρΔVDWWn]. (29)

Using the commutativity of and , we have and thus see that (29) is identical to (10).

### 4.2 No-Slip Boundary Case: Hypothetical Scheme

When Scheme 1 is applied to the no-slip boundary case with replaced by (see (9)), it does not exactly solve (10). Since is solved with the homogeneous Dirichlet boundary condition and is solved with the homogeneous Neumann boundary condition, the resulting does not necessarily satisfy the homogeneous Dirichlet boundary condition. While its normal component is zero, its parallel component is not guaranteed to be zero.

If we could impose the Dirichlet boundary condition to , the resulting would satisfy the homogeneous Dirichlet boundary condition. By expressing the vector Laplacian operator with a specified Dirichlet boundary condition in terms of and (see (9)), this procedure can be written as

 ~U=Un+νΔtL0(Un+~U2)+√2kBTνΔtρΔVDWWn+12νΔtBGΦ, (30a) DGΦ=D~U, (30b) Un+1=~U−GΦ. (30c)

However, the first two steps in this scheme cannot be sequentially implemented because is not available when is computed in (30a) and only available after (30b). Nonetheless, this hypothetical scheme is worth investigating because it exactly solves the Crank–Nicolson update (10).

We show that (30) is equivalent to (10) by using the commutativity of and , i.e.

 P(L0+B)=(L0+B)P. (31)

For the proof of (31), see B. By applying to (30a) and using and where , we obtain

 (32)

Since it can be shown using (31) that

 PL0~U=PPL0~U=P(L0P~U+BP~U−PB~U)=PL0Un+1−PBQ~U, (33)

(32) becomes identical to (10) and therefore the hypothetical scheme (30) would not introduce any splitting error.

### 4.3 No-Slip Boundary Case: K-Iteration Scheme (1≤K<∞)

To construct implementable schemes based on the hypothetical scheme (30), we first observe that satisfies the following fixed-point problem:

 ~U=(I−12νΔtL0)−1[(I+12νΔtL0)Un+12νΔtBQ~U+√2kBTνΔtρΔVDWWn]≡Ψ(~U). (34)

We then construct the following iteration procedure of computing , , , , to obtain the convergent solution and :

 ~Uk=Un+νΔtL0(Un+~Uk2)+√2kBTνΔtρΔVDWWn+12νΔtBGΦk−1, (35a) DGΦk=D~Uk, (35b)

where we assume . Since , it is easy to see that, if (35) converges, the limit is the solution of the fixed-point problem (34) and, equivalently, the solution of (30a) and (30b) in the hypothetical scheme. Therefore, the no-slip boundary solution for the next timestep is obtained as . We also note that the last term in (35a) is the boundary condition correction using .

Since the iteration (35) can be written as , where

 M=12νΔt(I−12νΔtL0)−1BQ, (36)

the convergence criterion is that the spectral radius of is smaller than 1. The rate of convergence (i.e. ) is given as

 r=ρ(M). (37)

Finally, by specifying a finite number of iterations (, we obtain Scheme 2, which we call the -iteration scheme. Note that the 1-iteration scheme corresponds to the projection method where no boundary correction is considered. Alternatively, one can impose a convergence criterion such as

 ∥Φk−1−Φk∥∥Φk∥≤ϵ. (38)

## 5 Steady-State Covariance Error Analysis

In this section, we analyze the accuracy of the -iteration scheme (see Scheme 2) by investigating the resulting steady-state covariance matrix . While the same results can be obtained using the structure factor (note that the structure factor is basically the Fourier transform of the steady-state covariance), we use the steady-state covariance in this section; we analyze the structure factor in section 6, where numerical validation results are presented.

As shown in section 4.2, if infinite iterations were performed each timestep, the resulting -iteration projection method would give the identical temporal update as the Crank–Nicolson scheme and thus its steady-state covariance would not have any temporal integration errors (i.e. ). When a finite number of iterations are used, the new state is computed from inexact , causing temporal integration errors in . The main theoretical result of this section is that the temporal integration error of the -iteration scheme in the steady-state covariance is of the order of :

 R(K)Δx,Δt−RΔx=O(ΔtK). (39)

In section 5.1, we show that temporal integration errors committed at each timestep due to a finite number of boundary corrections are :

 ~UK−~U∞=O(ΔtK). (40)

We note that this result strongly supports (39). In section 5.2, we confirm (39) using a semi-analytic approach. That is, by noting that can be determined as the unique solution of a linear system described in A, we directly compute it for specific values of and for some small-sized systems by solving the linear system.

### 5.1 Analytic Results

To show (40), we first derive expressions of and so that the temporal update of the -iteration scheme can be expressed in the compact form

 Un+1=P~UK=P[AKUn+BKWn]. (41)

By introducing , we express as

 ~U1=A−1−[A+Un+√2kBTνΔtρΔVDWWn]. (42)

Since , where , we recursively express for , and obtain the following general expression:

 ~UK=A−1−[K−1∑k=0(12νΔtBQA−1−)k][A+Un+√2kBTνΔtρΔVDWWn]. (43)

Using the identity , we have an alternative expression for , which gives the following expressions for and in the temporal update (41):

 AK=[K−1∑k=0(12νΔtA−1−BQ)k]A−1−A+,BK=√2kBTνΔtρΔV[K−1∑k=0(12νΔtA−1−BQ)k]A−1−DW. (44)

Using the well-known result for the geometric series of a matrix, we finally obtain

 ~UK=[I−(12νΔtA−1−BQ)K]~U∞ (45)

and thus (40).

We note that no notion of stochastic accuracy (e.g. weak vs. strong) is required to interpret (40) or (45) since and are fixed during the boundary correction iterations. However, to define the order of temporal integration errors committed at each timestep, a notion of stochastic order of convergence is needed. In this paper, instead of analyzing the weak or strong orders of accuracy of our schemes, we focus on the convergence of the resulting steady-state covariance (and equivalently, the structure factor). For discussion of stochastic accuracy of FHD schemes, we refer the reader to DelongGriffithVandenEijndenDonev2013 .

### 5.2 Semi-Analytic Results

As described in A, one can derive a linear system ((63) and (66), or equivalently, (70)) that uniquely determines the covariance matrix, using the definition of discretized operators and the expressions of and (see (41) and (44

)). However, the solution cannot be given explicitly. Hence, we construct the linear system for a spatially discretized system of a specific size and compute its solution numerically. Although the solution has numerical errors due to floating-point arithmetic, these errors can be controllable and kept small comparable to machine precision. In this sense, this approach gives semi-analytic results, which should not be confused with stochastic simulation results (given in section

6) containing sampling errors.

We present results obtained from a two-dimensional system with cells, where no-slip boundary conditions are imposed on all boundaries. We note that the choice of the number of cells is arbitrary, and the conclusion should not change for moderate to large system sizes (roughly speaking, ) as we justify below. However, we point out that semi-analytic results for a larger system quickly becomes computationally infeasible. This is because the size of the extended linear system (70) has a matrix with components for a system domain with cells. We assume and define the dimensionless number

 β=νΔtΔx2 (46)

and investigate how the errors change as the value of varies. For faces and , we define the error at the -component as

 Ef,f′=ΔV(R(K)Δx,Δt)f,f′−kBTρPf,f′, (47)

and define the maximum error as

 εmax=max(f,f′)|Ef,f′|. (48)

Since is linearly proportional to , we simply set .

Figure 3 shows how the errors in the steady-state covariance matrix are distributed for the 2-iteration scheme with . In panel (a), the error matrix is displayed as a grayscale image where the brightness of a pixel represents the magnitude of each component . The image shows that errors are dominant along the diagonal (i.e. for ). This is because diagonal components of the steady-state covariance matrix tend to be larger than off-diagonal components. In panel (b), the magnitude of each diagonal component is shown at the corresponding face in the physical domain of the system. The image indicates that errors are dominant near the boundaries. This observation is consistent with the fact that the inexact boundary correction causes temporal integration errors of the -iteration scheme due to a finite number of iterations. We note that the essentially same error pattern is observed for different system sizes. For cells with , the error distributions at the corners of the domain and at the sides of the domain remain the same. Moreover, compared with the maximum error for , the corresponding values for , 12, 13 have negligible deviations ( 0.2%). Hence, we expect our observations for the system to remain valid for larger systems.

In Figure 4, we show the dependence of the maximum error on the number of iterations and the stability condition number . Figure 4a demonstrates that the steady-state covariance matrix obtained by the -iteration scheme has the accuracy (see (39)). The semi-log plot of versus in Figure 4b indicates that for a given value of the error decreases exponentially with respect to . Moreover, we confirm that the errors asymptotically decay with the rate of convergence given in (37) for large values of .

## 6 Numerical Validations

In this section, we present numerical validation results of our -iteration scheme. We solve the two-dimensional stochastic Stokes equation (2) with no-slip boundary conditions using the -iteration scheme and calculate the equilibrium structure factors and using the time trajectories of . We compare these numerical results with the exact results for the discretized system as well as the semi-analytical results for (available only for small systems).

We recall that the equilibrium structure factors are defined as

 SUxΔx,Δt=Vlimn→∞⟨^Ux(nΔt)^U∗x(nΔt)⟩,SUyΔx,Δt=Vlimn→∞⟨^Uy(nΔt)^U∗y(nΔt)⟩, (49)

where is the volume of the system and is the Fourier transform of . As mentioned in section 3, due to the no-slip boundary conditions, we consider the extended domain (see Figure 11) to define the Fourier modes. As a result, for a system with cells, there are Fourier modes and the Fourier spectrum is symmetric with respect to and (see Figure 5 for the case).

Using a simulated time trajectory, the equilibrium structure factors are calculated as follows. Since the procedure is exactly the same for , we only explain the case of . When the time trajectory is computed up to timesteps, we compute the following time average to estimate the ensemble average , where the first timesteps are discarded to not include non-stationary data:

 SUxΔx,Δt≈VN2−N1N2∑n=N1+1^Ux(nΔt)^U∗x(nΔt). (50)

Hence, a sufficiently large value of should be used to reduce the systematic error, whereas the value of

should be large enough to control the level of the statistical error. By the central limit theorem, the magnitude of the statistical error is asymptotically proportional to

. The exact structure factor of the spatially discretized case can be computed using given in (18). Similarly, the theoretical values of the numerical structure factor , which one would obtain from (50) in the limits and can be computed using the semi-analytic result .

We present here the numerical structure factor results of that were calculated using the -iteration scheme (i.e. ) for three different system sizes, , , and cells. Parameter values, , , and , were used. The initial velocity field was set to zero and the trajectory was calculated up to timesteps for the and cases and for the case. To compute the equilibrium structure factor, the first timesteps were discarded.

Figure 5 shows the numerical structure factor for the smallest system size. Since the semi-analytic result of is available for this case, we compare the simulation error with the theoretically expected error (i.e. without statistical errors) for validation purposes. We see that the time integration error in the equlibrium structure factor due to inexact boundary corrections is reasonably small for rather large even when one boundary correction is used per timestep. In fact, the plot of is visually indistinguishable from that of the exact result (i.e. for ).

It is instructive to observe some features of . First, for the no-slip boundary and periodic boundary cases, their equilibrium structures are overall similar but different. For the periodic boundary case, is given as (see (16))

 SUxΔx=kBTρk2yk2x+k2y, (51)

and thus along a ray (i.e. ) the values of do not change. However, our no-slip boundary case result shows that the values of slightly change along a ray for larger values of and . Second, we see that becomes zero for the Fourier modes with . This is because if is independent of , it must be zero due to the boundary condition. In addition, also becomes zero for the Fourier modes with because these modes are omitted by the projection operator .

Figure 6 shows the results of the larger system sizes. As in the case, the plots of are visually indistinguishable from those of . As the number of cells increases, the plot of becomes similar to that of the continuum structure factor. The characteristic pattern observed in the error plot of the case appears in the error plots of the larger systems. Due to the smaller value of for the case, the level of statistical errors is relatively significant in the error plot. However, even for this case, the level of statistical errors in the structure factor plot is completely negligible.

## 7 Simulations of Giant Fluctuations

In this section, we apply our projection method to simulate the phenomenon of giant fluctuations. As experimentally observed in space GFexp ; GFexp2 , random advection (due to thermal fluctuations) can induce long-ranged concentration fluctuations when coupled with a concentration gradient in a micro-gravity environment. Specifically, in the absence of gravity, the nonequilibrium enhancement in the structure factor of the concentration fluctuations exhibits a power-law divergence, , where is the wavenumber. The theoretical explanation of the phenomenon OrtizDeZarateSengers2006 is regarded as one of the most significant accomplishments of the FHD approach. The incompressible computer simulation of the phenomenon was first performed in BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 , where the velocity equation was solved using the saddle-point system. In this section, we simulate giant fluctuations with similar settings considered in BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 . The goal of this section is to demonstrate that our projection method approach is readily applicable to the velocity equation coupled with the concentration equation.

### 7.1 Governing Equations

Following BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 , we consider a two-dimensional incompressible fluid system confined between two walls in the absence of gravity (see Figure 7). The fluid is a dilute solution, where the concentrations at the walls are held fixed with slightly different values. As a result the mean concentration profile has a small concentration gradient. The governing equations for the velocity and the mass fraction of the solute are written as

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩∂u∂t=νΔu−∇π+√2kBTνρ∇⋅Wv,∇⋅u=0,∂c∂t+u⋅∇c=χΔc+√2c(1−c)Mχρ∇⋅Wc, (52)

where is the constant density of the solution, is the kinematic viscosity of the solution, is the diffusivity of the solute in the solution, and is the mass of a solute molecule. and are independent spatiotemporal Gaussian white noise tensor fields. Note that we assume that viscous effects dominate inertial effects and omit the term. The size of the system domain is . At the walls, Dirichlet boundary conditions are imposed for : at , where is the mean mass fraction and . For the velocity, no-slip boundary conditions are imposed on the walls. Periodic boundary conditions are imposed for and in the horizontal direction.