Fluctuating hydrodynamic modeling of fluids at the nanoscale

A good representation of mesoscopic fluids is required to combine with molecular simulations at larger length and time scales [De Fabritiis et al., Phys. Rev. Lett. 97, 134501 (2006)]. However, accurate computational models of the hydrodynamics of nanoscale molecular assemblies are lacking, at least in part because of the stochastic character of the underlying fluctuating hydrodynamic equations. Here we derive a finite volume discretization of the compressible isothermal fluctuating hydrodynamic equations over a regular grid in the Eulerian reference system. We apply it to fluids such as argon at arbitrary densities and water under ambient conditions. To that end, molecular dynamics simulations are used to derive the required fluid properties. The equilibrium state of the model is shown to be thermodynamically consistent and correctly reproduces linear hydrodynamics including relaxation of sound and shear modes. We also consider nonequilibrium states involving diffusion and convection in cavities with no-slip boundary conditions.

Figure 1

(Color online) Velocity field $v_x$ as a function of $z$ (in Å) for this decaying transversal wave. Diamonds correspond to simulation results and the continuum lines correspond to the theoretical profile at snapshots corresponding to times $t=0,50,100,\dots ,750\mathrm{ps}$. At time $t=0\mathrm{ps}$ the amplitude is maximum while at $t=750\mathrm{ps}$ the amplitude is minimum.

Figure 2

(Color online) Relative error for the shear viscosity, sound absorption coefficient, and sound velocity as a function of the spatial resolution $\delta$. For any fluid parameter (e.g., $\nu$) the relative error is defined as $E_{\nu }\equiv \mid {\nu }_{\mathrm{num}}-\nu \mid ∕\nu$, where $\nu$ is the input value and ${\nu }_{\mathrm{num}}$ that measured from the relaxation of the corresponding hydrodynamic mode.

Figure 3

(Color online) The mass density $\rho$ as a function of $x$ (in Å) for a decaying longitudinal wave. Diamonds correspond to simulation results and the continuum lines correspond to the theoretical profile given by Eq. (12). The snapshots correspond to the times $t=0,25,50,\dots ,750\mathrm{ps}$. The agreement is remarkable. At time $t=0\mathrm{ps}$ the amplitude is maximum while at $t=750\mathrm{ps}$ the amplitude is minimum.

Figure 4

The equilibrium distribution of the $x$ component of the velocity at one fluid cell of volume $37.5{\mathrm{nm}}^3$ in a simulation of argon (circles) at $T=300\mathrm{K}$ and $\rho =0.6\mathrm{g}∕\mathrm{mol}∕{\mathrm{\AA }}^3$ compared with the theoretical normal distribution (continuum line). In this simulation $\mathrm{\Delta }t=20\mathrm{fs}$ and the best normal fit to the numerical distribution yields $T_{\mathrm{num}}=296.28\mathrm{K}$, that is a relative error of around 1.2%.

Figure 5

Relative error in the temperature of the scheme, $E_T=\mid T_{\mathrm{num}}-T\mid ∕T$, with $T_{\mathrm{num}}\equiv {\sum }_{\alpha }\mathrm{var}\left[v_{\alpha }\right]∕\left(3k_B\right)$ and $T=300\mathrm{K}$ the input equilibrium temperature. Results correspond to three-dimensional simulations of water at $\rho =0.6\mathrm{g}∕\mathrm{mol}∕{\mathrm{\AA }}^3$ with spatial resolution $\delta =20\mathrm{\AA }$ in each direction.

Figure 6

(Color online) Standard deviation of the mass density $S\left[\rho \right]$ in a fluid cell with volume $37.5{\mathrm{nm}}^3$ of argon and water at $T=300\mathrm{K}$. The continuous line corresponds to the grand canonical results obtained using the equation of state for argon, the dotted lines for water (TIP3P model) $[\rho k_{B}T∕(c_T^2{V}_c){]}^{1∕2}$, and the dashed lines show the ideal gas limits $\left[\right(m∕V)\rho {]}^{1∕2}$ with $m_{\mathrm{Ar}}=39.498\mathrm{g}∕\mathrm{mol}$ (upper curve) and $m_{{\mathrm{H}}_2\mathrm{O}}=18.015\mathrm{g}∕\mathrm{mol}$ (lower curve). Circles are results from the fluctuating hydrodynamics solver (using a water simulation at $\rho =0.632\mathrm{g}∕\mathrm{mol}∕{\mathrm{\AA }}^3$).

Figure 7

(Color online) Water at $\rho =0.632\mathrm{g}∕\mathrm{mol}∕{\mathrm{\AA }}^3$ and $T=300\mathrm{K}$ (ambient conditions) within a periodic box of size $50\times 50\times 2250{\mathrm{\AA }}^3$ at equilibrium; the FH mesh is comprised of $1\times 1\times 150$ cells. (a) The density field in the real space, $\rho (x,t_0)$. (b) The time dependence of the density at one cell $\rho (x_0,t)$. (c) The Fourier mode $\rho (k_1,t)$ associated with wave vector ${\mathit{k}}_1=(0,0,2\pi ∕L_z)$. (d) The (normalized) time correlation function $⟨\rho (k_1,t)\rho (k_1,0)⟩$. Note the different time scales associated with variations of quantities in (b), (c), and (d). The entire run is of $40\mathrm{ns}$ duration.

Figure 8

(Color online) The best fits to the decay times obtained from the time correlations of the Fourier modes of wavelength $\lambda =2\pi ∕k$. Symbols correspond to numerical results and lines to theoretical predictions [see Eqs. (14)]. (a) Shear decay times: the lines are the theoretical results $(k^2\nu {)}^{-1}$, (b) sound attenuation time $(k^2{\Gamma }_T{)}^{-1}$, and (c) the sound period $\lambda ∕c_T$. Circles are obtained from the imaginary part of $⟨h(k,\tau )h(k,0)⟩$ and squares from the real part, where $h=v_{\perp }$ (transversal velocity) in (a) while $h=\rho$ in (b) and (c) (similar results, not shown, were obtained for the longitudinal velocity $v_{\parallel }$). Results for water ($c_T=14.75\mathrm{\AA }∕\mathrm{ps}$, $\nu =84.98{\mathrm{\AA }}^2∕\mathrm{ps}$, and $\zeta ∕\rho =201.02{\mathrm{\AA }}^2∕\mathrm{ps}$) correspond to the simulations illustrated in Fig. 7, while argon $(c_T=5.61\mathrm{\AA }∕\mathrm{ps})$ corresponds to $\rho =0.6\mathrm{g}∕\mathrm{mol}∕{\mathrm{\AA }}^3$ and $T=300\mathrm{K}$ in the same box used in Fig. 7.

Figure 9

Boundary conditions imposed at the interface “$w$.” Cell 0 represents a fluid cell, and cells 1 and 2 are ghost cells. Arrows indicate the velocity vectors in a case with zero velocity imposed at “$w$” (no-slip boundary condition and rigid wall at rest).

Figure 10

(Color online) Stationary Couette profile according to fluctuating hydrodynamics: The diamonds are simulation results and the continuous line is the theoretical stationary linear profile. The wall velocity has been set to $2.04\mathrm{\AA }∕\mathrm{ps}$. The inset figure shows a deterministic Navier-Stokes simulation with the same parameters. Vertical dashed lines represent the boundary walls.

Figure 11

(Color online) Stationary Poiseuille profile according to fluctuating hydrodynamics: Diamonds correspond to simulation results and the continuous line is the theoretical profile associated with Navier-Stokes flow. The inset figure shows the same flow case for a purely deterministic simulation. In both figures, the vertical dashed lines represent the walls.

Figure 12

Streamline plots of a cavity flow for TIP3P water model in the stationary regime. The wall velocity along $y$ is $1\mathrm{\AA }∕\mathrm{ps}$, and the temperature is $300\mathrm{K}$. The dimensions of the cavity are $1500\times 1500\times 50{\mathrm{\AA }}^3$ and the mesh is $30\times 30\times 1$. In (a) thermal fluctuations of the pressure tensors are not present, and the Navier-Stokes solution recovered. In (b) we show the fluctuating hydrodynamics solution.

Figure 13

(Color online) Argon equation of state: Pressure (in bar) vs mass density at temperature $300\mathrm{K}$. The simulation results appear with error bars and the dashed line is the equation of state for the theoretical Lennard-Jones fluid in Refs. 30, 31, 32. The continuous line is the fit of the numerical data to the second order polynomial $3088.21-12065.2\rho +14765.8{\rho }^2$.

Figure 14

(Color online) Water pressure (in bar) vs temperature obtained from MD simulations of the TIP3P water model (using NAMD) at a fixed density $\rho =0.55066\mathrm{g}∕\left(\mathrm{mol}{\mathrm{\AA }}^3\right)$. The $*$ symbols correspond to simulations with all bonds rigid while the $◻$ symbols are for nonrigid bonds.

Figure 15

(Color online) Water pressure (in bars) versus mass density for a temperature of $300\mathrm{K}$ from MD simulations of the nonrigid TIP3P water model (results obtained with the NAMD code). The continuous line is the best fit of the numerical data to the second order polynomial $p\left(\rho \right)=38373.6-157398\rho +152881{\rho }^2$ bar.