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 kinesinmicrotubule systems lopes_membrane_2019 , nanofabricated devices for biomolecule applications B917759K , nano heat transport technologies nanoheatreview ; Wang2009 , brainonachip platforms 10.3389/fbioe.2019.00100 , and microsensing 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 subnanofluidic 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 PDEbased (as opposed to particlebased) approaches for homogeneous fluids. The Landau–Lifshitz Navier–Stokes equations (i.e. compressible FHD equations) were numerically solved for the onespecies CCSE_PRE2007 and binary mixture CCSE_ESAIM2010 cases as well as the multispecies case CCSE_PRE2014 , which was extended to include stochastic reactions CCSE_JCP2015 . The incompressible FHD formulation BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 was extended to quasiincompressible fluids by the low Mach number formulation for the binary mixture CCSE_CAMCOS2014 ; CCSE_CAMCOS2015 and muitispecies cases CCSE_PhysFluids2015 . The quasiincompressible 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 finitevolume 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 saddlepoint linear system, using the generalized minimal residual method (GMRES), an efficient variablecoefficient finitevolume Stokes solver was developed CaiNonakaBellGriffithDonev2014 . Several time integrators for FHD equations (e.g. semiimplicit 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 reactiondiffusion 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 EulerianLagrangian method Atzberger2011 ; wang_fluctuating_2018 , the boundary integral formulation BaoRachhKeavenyGreengardDonev2018 , and the immersed boundary approach CCSE_PRF2021 .
This paper presents a projectionbased 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 divergencefree 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 noslip 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 largescale simulations, we propose a projectionbased 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 steadystate 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 finitevolume 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:
(1a)  
(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 selfadvected 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 divergencefree velocity fields, we linearize equation (1) as
(2) 
Here, we assume that
is a spatiotemporal Gaussian white noise tensor field with independent components having covariance
(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 divergencefree, i.e. . This assumption implies that remains to be divergencefree for all .
In this paper, we consider the periodic boundary condition and the noslip 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 selfadjoint, 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(4) 
Therefore, it can be seen in (2) that the stochastic term is linked to the viscous term LandauLifshtz1959 by the fluctuationdissipation balance FluctuationDissipation .
2.2 Spatial Discretization
To spatially discretize the stochastic incompressible Stokes equation (2), we employ the staggeredgrid discretization on uniform mesh developed in BalboaUsabiagaBellDelgadoBuscalioniDonvFaiGriffithPeskin2012 . While we consider the twodimensional case in this paper for clarity, the threedimensional case is essentially the same. We denote the spatially discretized equation as
(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 noslip boundary condition.
2.2.1 Variables and Operators
According to the markerandcell 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 secondorder centered finite difference. At face , the Hodge gradient is also computed using the secondorder 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 secondorder 5point stencil.
The diagonal components of the stochastic stress tensor ( and ) are stored at the cell centers, while the offdiagonal 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:
(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(7) 
While
is simply the identity matrix in the periodic boundary case, we will see that
needs to be modified for the noslip 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 fluctuationdissipation 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 fluctuationdissipation 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 secondorder centered finite differences, the discrete stochastic divergence is constructed. Then it can easily be seen that the following discrete fluctuationdissipation balance
(8) 
is satisfied. The properties of the discretized operators and are important for computing the steadystate covariance in sections 3 and 5.
2.2.2 Boundary Conditions
In the presence of nonperiodic boundaries, the discrete operators defined above need to be modified near the boundaries. For the noslip 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 5point stencil, and can be expressed as
(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 fluctuationdissipation 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
(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 noslip boundary to satisfy the discrete fluctuationdissipation 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 :
(11a)  
(11b) 
The main goal of this paper is to solve (10) not using a linear solver for the saddlepoint system (11) but using a projection method.
3 SteadyState 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 steadystate 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 steadystate covariances are defined asContinuum  (12a)  
Spatially discretized  (12b)  
Fully discretized  (12c) 
For each case, the static structure factor is defined as the Fourier transform of the steadystate covariance. We note that, for the noslip 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 projectionmethodbased schemes by computing the steadystate 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 projectionmethodbased schemes. In A, we develop a systematic procedure to construct and solve linear systems from which the steadystate covariance can be uniquely determined.
3.1 Continuum and Spatially Discretized Cases
Using the result of A, the steadystate covariance of the continuum equation (2) is given as the unique solution of
(13a)  
(13b) 
The physical intuition that the equilibrium covariance must be proportional to the projection operator, suggests that
(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 steadystate covariance (with normalization DonevVandenEijndenGarciaBell2010 with respect to the volume of the system ):
(15) 
where a hat denotes the Fourier transform. For the periodic boundary case, the structure factor has a compact form
(16) 
It is easy to see from (16) that all divergencefree 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 fluctuationdissipation balance is given in the form of (8), the linear system that uniquely determines is
(17a)  
(17b) 
and, similarly to the continuum steadystate expression (14), the discrete steadystate covariance is given as
(18) 
3.2 Crank–Nicolson Scheme
Introducing , the Crank–Nicolson scheme (10) can be written as
(19) 
Hence, using the result of A, the steadystate covariance is given as the unique solution of
(20a)  
(20b) 
By observing that
(21) 
and using the discrete fluctuationdissipation balance (8), one can show that
(22) 
satisfies (20). The results (18) and (22) show that the Crank–Nicolson scheme does not introduce any temporal integration errors to the steadystate 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 noslip boundary cases are given in Schemes 1 and 2, respectively. Our projectionmethodbased 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 noslip case and introduce the idea of iterative boundary corrections. In section 4.3, we present our iterative scheme for the noslip case.
Given velocity at timestep , the velocity at the next timestep is updated via as follows.
(23) 
(24) 
(25) 
Given velocity at timestep , the velocity at the next timestep is updated as follows.

For :

Compute using . For , use .
(26) 
Compute .
(27)


Compute from the projection of :
(28)
4.1 Periodic Boundary Case
4.2 NoSlip Boundary Case: Hypothetical Scheme
When Scheme 1 is applied to the noslip 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
(30a)  
(30b)  
(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.
(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
(33) 
(32) becomes identical to (10) and therefore the hypothetical scheme (30) would not introduce any splitting error.
4.3 NoSlip Boundary Case: Iteration Scheme
To construct implementable schemes based on the hypothetical scheme (30), we first observe that satisfies the following fixedpoint problem:
(34) 
We then construct the following iteration procedure of computing , , , , to obtain the convergent solution and :
(35a)  
(35b) 
where we assume . Since , it is easy to see that, if (35) converges, the limit is the solution of the fixedpoint problem (34) and, equivalently, the solution of (30a) and (30b) in the hypothetical scheme. Therefore, the noslip 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
(36) 
the convergence criterion is that the spectral radius of is smaller than 1. The rate of convergence (i.e. ) is given as
(37) 
Finally, by specifying a finite number of iterations (, we obtain Scheme 2, which we call the iteration scheme. Note that the 1iteration scheme corresponds to the projection method where no boundary correction is considered. Alternatively, one can impose a convergence criterion such as
(38) 
5 SteadyState Covariance Error Analysis
In this section, we analyze the accuracy of the iteration scheme (see Scheme 2) by investigating the resulting steadystate 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 steadystate covariance), we use the steadystate 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 steadystate 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 steadystate covariance is of the order of :
(39) 
In section 5.1, we show that temporal integration errors committed at each timestep due to a finite number of boundary corrections are :
(40) 
We note that this result strongly supports (39). In section 5.2, we confirm (39) using a semianalytic 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 smallsized 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
(41) 
By introducing , we express as
(42) 
Since , where , we recursively express for , and obtain the following general expression:
(43) 
Using the identity , we have an alternative expression for , which gives the following expressions for and in the temporal update (41):
(44) 
Using the wellknown result for the geometric series of a matrix, we finally obtain
(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 steadystate covariance (and equivalently, the structure factor). For discussion of stochastic accuracy of FHD schemes, we refer the reader to DelongGriffithVandenEijndenDonev2013 .
5.2 SemiAnalytic 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 floatingpoint arithmetic, these errors can be controllable and kept small comparable to machine precision. In this sense, this approach gives semianalytic results, which should not be confused with stochastic simulation results (given in section
6) containing sampling errors.We present results obtained from a twodimensional system with cells, where noslip 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 semianalytic 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
(46) 
and investigate how the errors change as the value of varies. For faces and , we define the error at the component as
(47) 
and define the maximum error as
(48) 
Since is linearly proportional to , we simply set .
Figure 3 shows how the errors in the steadystate covariance matrix are distributed for the 2iteration 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 steadystate covariance matrix tend to be larger than offdiagonal 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.
is estimated from the spectral radius of the iteration matrix, see (
37).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 steadystate covariance matrix obtained by the iteration scheme has the accuracy (see (39)). The semilog 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 twodimensional stochastic Stokes equation (2) with noslip 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 semianalytical results for (available only for small systems).
We recall that the equilibrium structure factors are defined as
(49) 
where is the volume of the system and is the Fourier transform of . As mentioned in section 3, due to the noslip 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 nonstationary data:
(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 semianalytic 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 semianalytic 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 noslip boundary and periodic boundary cases, their equilibrium structures are overall similar but different. For the periodic boundary case, is given as (see (16))
(51) 
and thus along a ray (i.e. ) the values of do not change. However, our noslip 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 longranged concentration fluctuations when coupled with a concentration gradient in a microgravity environment. Specifically, in the absence of gravity, the nonequilibrium enhancement in the structure factor of the concentration fluctuations exhibits a powerlaw 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 saddlepoint 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 twodimensional 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
(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, noslip boundary conditions are imposed on the walls. Periodic boundary conditions are imposed for and in the horizontal direction.
Comments
There are no comments yet.