On the Suppression and Distortion of Non-Equilibrium Fluctuations by Transpiration

TL;DR

This paper investigates how a transpiration interface influences long-range hydrodynamic fluctuations in a non-equilibrium dilute gas, using both DSMC and fluctuating hydrodynamics simulations, revealing significant effects and computational advantages.

Contribution

It introduces models for the interface in fluctuating hydrodynamics and compares them with DSMC, demonstrating the impact on non-equilibrium fluctuation correlations.

Findings

Good agreement between DSMC and FHD results.

Transpiration significantly alters long-range correlations.

FHD offers computational speed advantages.

Abstract

A fluid in a non-equilibrium state exhibits long-ranged correlations of its hydrodynamic fluctuations. In this article, we examine the effect of a transpiration interface on these correlations -- specifically, we consider a dilute gas in a domain bisected by the interface. The system is held in a non-equilibrium steady state by using isothermal walls to impose a temperature gradient. The gas is simulated using both direct simulation Monte Carlo (DSMC) and fluctuating hydrodynamics (FHD). For the FHD simulations two models are developed for the interface based on master equation and Langevin approaches. For appropriate simulation parameters, good agreement is observed between DSMC and FHD results with the latter showing a significant advantage in computational speed. For each approach we quantify the effects of transpiration on long-ranged correlations in the hydrodynamic variables.

Tables

Table 1. Table 1 : Parameters used when comparing FHD and DSMC simulations, corresponding to a dilute argon gas, in cgs units. For DSMC simulations the molecular properties are applied directly, for FHD simulations they are used to calculate the viscosity and thermal conductivity (see Ref.  11 )

Molecular diameter	$3.66\times {10}^{-8}$		System length ($x$ dir)	$2.5\times {10}^{-4}$
Molecular mass	$6.63\times {10}^{-23}$		Interface $x$ location	$1.25\times {10}^{-4}$
Reference temperature	$273$		Cell size $\Delta x$	$3.125\times {10}^{-6}$
Sound speed	$30,781$		Time step $\Delta t$	$1.0\times {10}^{-12}$
Specific heat $c_v$	$3.12\times {10}^6$		Number of cells ($x$ dir)	80
Reference mean free path	$6.26\times {10}^{-6}$		System height, width ($y,z$ dir)	$1.25\times {10}^{-6}$
Reference mean free time	$1.64\times {10}^{-10}$		Average particles per cell	$131$

Table 2. Table 2 : Boundary temperatures extrapolated from DSMC results, for use in FHD simulations. The DSMC wall temperatures were 273 K and 519 K.

$f$		$T\left(0\right)$		$T\left(L\right)$
no membrane		280.8		508.5
0.1		277.6		509.9
0.05		276.2		515.4
0.025		275.0		516.9
0.0125		274.1		517.9
0.00625		273.6		518.4

Equations

$$\frac{\partial }{\partial t}\left(\rho \right)=-\nabla \cdot \left(\rho \mathbf{u}\right),$$

$$\frac{\partial }{\partial t}\left(\rho \mathbf{u}\right)=-\nabla \cdot \left(\rho \mathbf{u}\otimes \mathbf{u}\right)-\nabla P-\nabla \cdot \left[\mathbf{\Pi }+\tilde{\mathbf{\Pi }}\right],$$

$$\frac{\partial }{\partial t}\left(\rho E\right)=-\nabla \cdot \left(\rho \mathbf{u}E+p\mathbf{u}\right)-\nabla \cdot \left[\mathbf{Q}+\tilde{\mathbf{Q}}\right]-\nabla \cdot \left(\left[\mathbf{\Pi }+\tilde{\mathbf{\Pi }}\right]\cdot \mathbf{u}\right).$$

$${\Pi }_{ij}=-\eta \left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)+{\delta }_{ij}\left(\frac{2}{3}\eta \nabla \cdot \mathbf{u}\right),\mathrm{and}\mathbf{Q}=-\kappa \nabla T,$$

$$⟨{\tilde{\Pi }}_{ij}(\mathbf{r},t),{\tilde{\Pi }}_{kl}({\mathbf{r}}^{\prime },t^{\prime })⟩=2k_B\eta T\left[({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk})-\frac{1}{3}{\delta }_{ij}{\delta }_{kl}\right]\delta (t-t^{\prime })\delta (\mathbf{r}-{\mathbf{r}}^{\prime }),$$

$$⟨{\tilde{Q}}_i(\mathbf{r},t),{\tilde{Q}}_j({\mathbf{r}}^{\prime },t^{\prime })⟩=2k_B\kappa T^2{\delta }_{ij}\delta (t-t^{\prime })\delta (\mathbf{r}-{\mathbf{r}}^{\prime }).$$

$${\tau }_{\mathrm{e}}=\frac{m}{\rho fA}\sqrt{\frac{2\pi m}{k_{B}T}}.$$

$${\tau }_{\mathrm{c}}=\frac{m}{4\sqrt{\pi }\rho d^2}\sqrt{\frac{2\pi m}{k_{B}T}},$$

$$P_{ϵ}(ϵ)=\frac{ϵ}{k_B^2{T}^2}{e}^{-ϵ/(k_{B}T)};$$

$$W_{\to }(ϵ)=\frac{fA{n}_L}{\sqrt{2\pi m{k}_B{T}_L}}\frac{ϵ}{k_B{T}_L}{e}^{-ϵ/(k_B{T}_L)},$$

$$W_{\leftarrow }(ϵ)=\frac{fA{n}_R}{\sqrt{2\pi m{k}_B{T}_R}}\frac{ϵ}{k_B{T}_R}{e}^{-ϵ/(k_B{T}_R)},$$

$$\frac{d}{dt}P(M,\mathcal{E})={\int }_0^{\infty }P(M-m,\mathcal{E}-ϵ)$$

$$-$$

$$⟨J_M{⟩}_{\to }={\rho }_{L}fA\sqrt{\frac{k_B{T}_L}{2\pi m}},$$

$$⟨J_M⟩=⟨J_M{⟩}_{\to }-⟨J_M{⟩}_{\leftarrow }=\frac{fA{k}_B^{1/2}}{\sqrt{2\pi m}}\left({\rho }_L{T}_L^{1/2}-{\rho }_R{T}_R^{1/2}\right).$$

$$⟨J_{\mathcal{E}}⟩=\frac{2k_B}{m}\left(T_L⟨J_M{⟩}_{\to }-T_R⟨J_M{⟩}_{\leftarrow }\right)=\frac{2fA{k}_B^{3/2}}{m\sqrt{2\pi m}}\left({\rho }_L{T}_L^{3/2}-{\rho }_R{T}_R^{3/2}\right).$$

$$⟨J_{\delta M\delta M}⟩=\frac{mfA{k}_B^{1/2}}{\sqrt{2\pi m}}\left({\rho }_R{T}_R^{1/2}+{\rho }_L{T}_L^{1/2}\right),$$

$$⟨J_{\delta \mathcal{E}\delta \mathcal{E}}⟩=\frac{6fA{k}_B^{5/2}}{m\sqrt{2\pi m}}\left({\rho }_R{T}_R^{5/2}+{\rho }_L{T}_L^{5/2}\right),$$

$$⟨J_{\delta \mathcal{E}\delta M}⟩=\frac{2fA{k}_B^{3/2}}{\sqrt{2\pi m}}\left({\rho }_R{T}_R^{3/2}+{\rho }_L{T}_L^{3/2}\right).$$

$$\frac{d}{dt}{M}_R=⟨J_M⟩+{\tilde{J}}_M,$$

$$\frac{d}{dt}{\mathcal{E}}_R=⟨J_{\mathcal{E}}⟩+{\tilde{J}}_{\mathcal{E}},$$

$$\frac{d}{dt}{M}_L=-\frac{d}{dt}{M}_R,\frac{d}{dt}{\mathcal{E}}_L=-\frac{d}{dt}{\mathcal{E}}_R,$$

$$⟨{\tilde{J}}_M(t){\tilde{J}}_M(t^{\prime })⟩=⟨J_{\delta M\delta M}⟩\delta (t-t^{\prime }),$$

$$⟨{\tilde{J}}_{\mathcal{E}}(t){\tilde{J}}_{\mathcal{E}}(t^{\prime })⟩=⟨J_{\delta \mathcal{E}\delta \mathcal{E}}⟩\delta (t-t^{\prime }),$$

$$⟨{\tilde{J}}_{\mathcal{E}}(t){\tilde{J}}_M(t^{\prime })⟩=⟨J_{\delta \mathcal{E}\delta M}⟩\delta (t-t^{\prime }).$$

$$\theta =-\frac{\log ({\mathcal{R}}_1)}{r_L+r_R},$$

$$ϵ=-k_B{T}_{\ast }\log ({\mathcal{R}}_3{\mathcal{R}}_4),$$

$$r=\frac{⟨J_{\delta \mathcal{E}\delta M}⟩}{\sqrt{⟨J_{\delta \mathcal{E}\delta \mathcal{E}}⟩⟨J_{\delta M\delta M}⟩}}.$$

$${\mathcal{G}}_3={\mathcal{G}}_{1}r+\sqrt{1-r^2}{\mathcal{G}}_2$$

$$F_M=\Delta t⟨J_M⟩+\sqrt{\Delta t⟨J_{\delta M\delta M}⟩}{\mathcal{G}}_1,$$

$$F_{\mathcal{E}}=\Delta t⟨J_{\mathcal{E}}⟩+\sqrt{\Delta t⟨J_{\delta \mathcal{E}\delta \mathcal{E}}⟩}{\mathcal{G}}_3.$$

Full text

I Introduction

Recent research in the field of fluid dynamics has seen an increased focus on nano and microscale applications, and in particular, the role of thermal fluctuations. These fluctuations are important in areas ranging from microfluidic engineeeringHänggi and Marchesoni (2009) to molecular biologyMikhailov and Kapral (2015). Fluctuations also yield information about the transport properties of the fluid itselfBerne and Pecora (2013), and are important in dynamic processes such as phase transitions, combustion and ignition, nucleation, hydrodynamic instabilities, and a host of Brownian dynamics phenomena.Lemarchand and Nowakowski (2004); Donev et al. (2011a); Vailati et al. (2011); Morgado, Nowakowski, and Lemarchand (2018); Gallis et al. (2015); Kirkpatrick, Ortiz de Zárate, and Sengers (2015)

Traditional computational fluid dynamic approaches have focused on solutions to the deterministic Navier-Stokes equations, which neglect thermal fluctuations. To capture this important feature of nanoscale flows particle based methods are typically employed. These techniques operate by replicating, with varying degrees of approximation, the dynamics of the fluid’s constituent molecules. Molecular dynamicsFrenkel and Smit (2001) takes the approach of deterministically simulating the interactions between molecules, often using a relatively complex interaction potential. Although this approach accurately simulates fluids at the nanoscale, its computational cost makes it inappropriate for mesoscale applications where thermal fluctuations are of importance. The direct simulation Monte CarloBird (1994, 2013) (DSMC) method greatly simplifies the particle interactions using a stochastic approach, allowing the simulation of much larger systems while still reproducing the correct thermal fluctuations. DSMC solves the Boltzmann equationWagner (1992) so it can be used for flows of arbitrary Knudsen number,Sone (2000) whereas solutions of the Navier-Stokes equations are valid only when the Knudsen number is sufficiently small. However this technique applies only to dilute gasses, and is still relatively computationally expensive.

An alternative method for simulating mesoscopic systems is given by fluctuating hydrodynamics (FHD). Following the approach of Landau and Lifshitz Landau and Lifshitz (1959); de Zarate and Sengers (2006), FHD extends the deterministic Navier-Stokes equations by incorporating stochastic fluxes that lead to fluctuations that are consistent with statistical mechanics. The validity of fluctuating hydrodynamics for non-equilibrium systems has been confirmed experimentally by light scatteringSegrè et al. (1992) and shadowgraph experimentsTakacs et al. (2011). Numerical methods for solving the compressible and incompressible FHD equations are discussed in Refs. 19; 20; 21, and the approach has been further extended to multispecies mixtures,Bell, Garcia, and Williams (2010); Donev et al. (2011a, b); Balakrishnan et al. (2014); Donev et al. (2014); Nonaka et al. (2015); Donev et al. (2015) multiphase flows,Chaudhri et al. (2014) reacting flows,Bhattacharjee et al. (2015) and electrokinetic flows.Péraud et al. (2016, 2017). Hybrid methods combining FHD and DSMC are investigated in Refs. 32; 33.

When a fluid is held in a non-equilibrium state, long range correlations in the thermal fluctuationsde Zarate and Sengers (2006) can result in macroscopic effects.Donev et al. (2014) These non-equilibrium fluctuations have been studied in many scenarios in which the correlations extend through the bulk of the fluid, however less is known regarding the effect of an interface, such as a porous wall or a thin membrane. This type of interface has previously been studied using molecular dynamicsCleuren, Van den Broeck, and Kawai (2006) in the context of two separate reservoirs of gas, each held in equilibrium at different temperatures and densities. In this article we consider a dilute gas in a volume bisected by simple porous transpiration interface, with the system driven out of equilibrium by inducing a temperature gradient. The partitioning interface is taken to be a thin adiabatic wall, with holes that are small compared with the mean free path of the gas. This ensures that the pores themselves do not have a significant hydrodynamic effect, yielding a simple system for the study of correlations. Particles passing through the interface conserve their energy, but not their momentum, consistent with the particles experiencing randomizing elastic collisions as they pass through the membrane.

In this article, both DSMC and FHD are used to simulate this system, yielding information on long range correlations and validating the use of FHD for transpiration interfaces. In Section II we begin with a brief review of DSMC, and of the theory underlying fluctuating hydrodynamics including the numerical implementation of the FHD equations.In Section III we present several models for transpiration interfaces, using both master equation and Langevin approaches in the FHD context, and the kinetic model applicable to DSMC simulations. In Section IV, we validate these approaches by comparing the long range correlations occurring in both FHD and DSMC simulations, in the presence of interfaces with varying porosity. The role of the interface in distorting and suppressing these correlations is discussed, and a brief comparison of the relative computational performance of the DSMC and FHD simulations is made. Some concluding remarks are given in Section V.

II Fluid Simulation Models

For the study of non-equilibrium fluctuations we use two independent approaches for the bulk fluid: a particle model, Direct Simulation Monte Carlo (DSMC), and a continuum model, fluctuating hydrodynamics (FHD). This section summarizes these numerical models; the treatment of the transpiration interface is described in Section III.

II.1 Direct Simulation Monte Carlo (DSMC)

Direct Simulation Monte Carlo is a well-known method developed by Graeme Bird for computing gas dynamics from the perspective of kinetic theory.Bird (1994) As in molecular dynamics, the state of the system in DSMC is given by the positions and velocities of particles. At each time step the particles are first independently moved ballistically, imposing boundary conditions if they strike a surface. Generally, isothermal boundary conditions are implemented using a diffuse surface interaction while adiabatic boundaries are simulated using specular reflection; surfaces with partial accommodation can also be simulated. Collisions are evaluated by a stochastic process, selecting the post-collision angles from their kinetic theory distributions while conserving momentum and energy. We employ the hard sphere collision operator in the present work but many other collisions models are available in DSMC. Pedagogical explanations of DSMC are contained in Refs. 35, 36, and 37 with a complete description given in Ref. 12.

Although DSMC is a stochastic algorithm the “Monte Carlo” elements of the method have nothing to do with the statistical variation of the hydrodynamic variables; the same fluctuations also occur in deterministic methods, such as molecular dynamics. Under a relatively weak set of conditions (e.g., detailed balance) discrete event processes reproduce the correct physical fluctuations at hydrodynamic scales.Keizer (1978, 2012) When each DSMC particle represents one molecule in the gas, the magnitude of the hydrodynamic fluctuations is in agreement with equilibrium statistical mechanics.Hadjiconstantinou et al. (2003) For both equilibrium and non-equilibrium systems DSMC yields the physical spectra of spontaneous fluctuations as predicted by theoryMansour et al. (1987); Garcia and Penland (1991) and observed in experimentsBruno et al. (2006); Bruno, Frezzotti, and Ghiroldi (2017). In this paper the numerical measurements from DSMC simulations are used to validate the fluctuating hydrodynamic simulation results.

II.2 Fluctuating Hydrodynamics (FHD)

The equations of fluctuating hydrodynamics are obtained by adding white noise to each dissipative flux in the deterministic Navier-Stokes equations,Giovangigli (2012) to represent thermal fluctuations. We consider the compressible fluctuating hydrodynamic equations for a pure monatomic ideal gas; generalizations are found in Refs. 22; 5; 23; 24; 25; 26; 27; 28; 29; 30.The continuity, momentum and energy conservation equations for FHD are

[TABLE]

[TABLE]

[TABLE]

Here, $\rho$ is the mass density of the gas, $\bf{u}$ is the velocity, $P$ is the pressure, $T$ is temperature, and $E=c_{V}T+\frac{1}{2}|{\bf u}|^{2}$ is the total specific energy. For a monatomic ideal gas the specific heat at constant volume is $c_{V}=\frac{3}{2}k_{B}/m$, where $k_{B}$ is Boltzmann’s constant and $m$ is the molecular mass; the equation of state is $P=nk_{B}T$ where $n=\rho/m$ is the number density. From Eqs. (1)-(3) we define the fluxes in the mass, momentum, and total energy as $J_{M}$, $J_{\mathbfcal{P}}$, and $J_{\mathcal{E}}$, respectively.

In the momentum and energy equations (Eqs. (2) and  (3)), the deterministic viscous stress tensor ${\bf\Pi}$ and deterministic heat flux vector ${\bf Q}$ are given by

[TABLE]

where $\delta_{ij}$ is the Kronecker delta, $\eta(T)$ is the shear viscosity, and $\kappa(T)$ is the thermal conductivity; consistent with a monatomic gas, we take the bulk viscosity to be zero.

The terms that appear with a tilde in Eqs. (2) and (3) denote stochastic contributions from Gaussian random fields. The stochastic viscous tensor, $\mathbf{\widetilde{\Pi}}$, stochastic heat flux tensor, $\mathbf{\widetilde{Q}}$, have zero mean and covariances Landau and Lifshitz (1959); de Zarate and Sengers (2006)

[TABLE]

and

[TABLE]

where $\delta(..)$ is the Dirac delta function.

The FHD Eqs. (1)-(3) above, are simulated numerically in one dimension using a finite volume method with third-order Runge-Kutta time stepping. A complete description of this approach is given in Sec VI of Ref. 20.

III Transpiration Interface Models

As discussed in Section I, we consider an adiabatic, stationary interface that is porous to the gas molecules (e.g., a specular wall with small holes). Particles encountering the interface are allowed to pass with an effusion probability, $f$. This probability can be viewed as the fraction of the wall area comprised of holes. The conditions in the gas are chosen such that a particle that crosses the interface experiences a large number of inter-molecular collisions before crossing again – this is called the effusion conditionBird (1994); Liepmann (1961).

III.1 DSMC interface model

The implementation of the transpiration interface is straightforward in DSMC. Particles reaching the interface are transmitted with probability $f$ and specularly reflected with probability $1-f$. The transmitted particles are assumed to interact with the adiabatic porous interface such that the particle’s energy, $\epsilon$, is conserved but not its momentum. Therefore the velocity of a transmitted particle, $\mathbf{v}$, is reset by randomly selecting it from a hemisphere of radius $|\mathbf{v}|=\sqrt{2\epsilon/m}$.

III.2 FHD interface models

The specular condition for the surface of the interface, discussed above, is modeled hydrodynamically as a full slip adiabatic boundary condition. We present two methods for modeling the effusion of mass and energy through the interface, based on master equation and Langevin approaches; as discussed above, the effusive momentum flux is set to zero. We first discuss the kinetic theory underlying both these techniques.

III.2.1 Kinetic theory

For a gas at equilibrium, the time between particle crossings through a transpiration interface of area $A$ is exponentially distributed, with the mean time between crossings given byBird (1994)

[TABLE]

By comparison, the mean time between collisions for a hard sphere particleBird (1994) is

[TABLE]

where $d$ is the collisional diameter of the particle. The effusion condition described above is $nA\lambda_{\mathrm{c}}\tau_{\mathrm{e}}\gg\tau_{\mathrm{c}}$ where $\lambda_{\mathrm{c}}$ is the mean free path, which is satisfied when $f\ll 1$. The probability distribution for the kinetic energy, $\epsilon$, of particles passing through the interface is,

[TABLE]

this is a gamma distribution, $\Gamma(2,k_{B}T)$.Gentle (2013) Eqs. (7) and (9) form the basis of the stochastic simulation algorithm given in Section III.2.2.

From the above we find that the probability per unit time of observing a particle of energy $\epsilon$ crossing the interface is

[TABLE]

[TABLE]

where the arrows indicate the direction of the crossing, and the subscripts $L$ and $R$ refer to the gas on the left and right hand sides of the interface. This leads to the master equation for $P(M,\mathcal{E})$, i.e., the probability density of the state on the right of the interface with mass $M$ and total energy $\mathcal{E}$,

[TABLE]

Following from the above, Cleuren, Van den Broeck, and Kawai (2006) the average one-way flux of mass passing from left to right through the interface is,

[TABLE]

with the arrow indicating the direction of flux. The total average mass flux is given by

[TABLE]

The average energy of a particle crossing the interface is $2k_{B}T$, so the average total energy flux is

[TABLE]

As discussed above, we assume that there is no momentum flux across the interface so $J{P}=0$. Since the non-equilibrium conditions considered in this article result in zero average net momentum transport, this difference in treatment between DSMC and FHD appears to have a negligible impact on the long-ranged correlations of interest (see Section IV.3).

The variances and covariance of the mass and energy fluxes are

[TABLE]

[TABLE]

[TABLE]

The resulting Langevin equations are given by

[TABLE]

[TABLE]

[TABLE]

where $\widetilde{J}_{M}$ and $\widetilde{J}_{\mathcal{E}}$ are Gaussian noise terms, with

[TABLE]

[TABLE]

[TABLE]

These equations form the basis of the Langevin method discussed in Section III.2.3.

III.2.2 Master equation implementation

Following from the theory presented above, we define a stochastic simulation algorithm (SSA)Gillespie (2007) implementation of the interface. This operates by repeatedly generating particle crossing times from an exponential distribution, until the cumulative sum of these crossing times exceeds a single FHD time step, $\Delta t$. The energy of each crossing particle is selected from the appropriate gamma distribution, and a tally of the total mass and energy transfer is used to update the hydrodynamic cells on each side of the interface.

The implementation of the master equation model over a single time step is as follows:

• •

Take $\rho_{L}$ and $\rho_{R}$ to be the density at the cell centers to the left and to the right of the interface. Similarly, the temperature and specific energy in said cells is indicated by $T_{L}$, $T_{R}$, $E_{L}$, and $E_{R}$.

• •

From Eq. (7) the rate of particles crossing from each direction is $r_{L}=1/\tau_{\mathrm{e}}(\rho_{L},T_{L})$ and $r_{R}=1/\tau_{\mathrm{e}}(\rho_{R},T_{R})$.

• •

Initialize the accumulation of mass and energy flux, $F_{M}$ and $F_{\mathcal{E}}$, across the interface to zero.

• •

Set the interface time counter $t_{\mathrm{ssa}}=0$. Loop while $t_{\mathrm{ssa}}<\Delta t$:

• –

Given an exponential distribution of time between particle crossings, the time until the next particle crossing is given by

[TABLE]

where $\mathcal{R}_{i}\in(0,1]$ indicates an independent uniform random number.

• –

Increment the interface time counter, $t_{\mathrm{ssa}}\mathrel{+}=\theta$. If $t_{\mathrm{ssa}}\geq\Delta t$ then break out of the loop.

• –

Select the direction for the particle crossing the interface. If $r_{L}/(r_{L}+r_{R})>\mathcal{R}_{2}$, then the particle crosses from left to right, otherwise it crosses from right to left.

• –

Following from Eq. (9), the energy of the effusing particle is gamma distributed and generated byGentle (2013)

[TABLE]

where $T_{*}=T_{L}$ if the particle is traveling left to right, or $T_{*}=T_{R}$ otherwise.

• –

If the particle is traveling left to right, $F_{M}\mathrel{+}=m$ and $F_{\mathcal{E}}\mathrel{+}=\epsilon$. Otherwise $F_{M}\mathrel{-}=m$ and $F_{\mathcal{E}}\mathrel{-}=\epsilon$.

• •

Finally, the left cell is updated as $\rho_{L}\mathrel{-}=F_{M}/\Delta V$ and $\rho_{L}E_{L}\mathrel{-}=F_{\mathcal{E}}/\Delta V$, and the right cell is updated as $\rho_{R}\mathrel{+}=F_{M}/\Delta V$ and $\rho_{R}E_{R}\mathrel{+}=F_{\mathcal{E}}/\Delta V$ with $\Delta V$ being the cell volume.

Note that the interface is modeled as a slip wall since the flux of momentum crossing the interface is zero.

III.2.3 Langevin equation implementation

In contrast to the SSA approach, which considers discrete particle crossing events, the Langevin algorithm generates the total mass and energy transfer based on the expected values. A Gaussian approximation is used for the fluctuations in these quantities – this is similar in principle to the stochastic fluxes generated between cells in the FHD algorithm. The implementation of the Langevin model over a single time step is as follows:

• •

Take $\rho_{L}$ and $\rho_{R}$ to be the density at the cell centers to the left and to the right of the interface. Similarly, the temperature and specific energy in said cells is indicated by $T_{L}$, $T_{R}$, $E_{L}$, and $E_{R}$.

• •

Calculate the mean mass and energy fluxes given by Eqs. (14) and (15), and the variances and covariance given by Eqs. (16)-(18).

• •

Calculate the correlation coefficient, given by

[TABLE]

• •

Generate correlated random variables,

[TABLE]

where $\mathcal{G}_{1}$ and $\mathcal{G}_{2}$ are independent Gaussian random numbers of mean zero and variance one.

• •

Following from Eq. (19), the mass and energy flux accumulations are given by

[TABLE]

[TABLE]

• •

Finally, the left cell is updated as $\rho_{L}\mathrel{-}=F_{M}/\Delta V$ and $\rho_{L}E_{L}\mathrel{-}=F_{\mathcal{E}}/\Delta V$, and the right cell is updated as $\rho_{R}\mathrel{+}=F_{M}/\Delta V$ and $\rho_{R}E_{R}\mathrel{+}=F_{\mathcal{E}}/\Delta V$ with $\Delta V$ being the cell volume.

Again, note that the interface is modeled as a slip wall since the flux of momentum crossing the interface is zero.

IV Simulation Results

In this section, we begin by validating the use of FHD for transpiration by comparison of the FHD and DSMC simulations. For the majority of this section, the master equation method of Section III.2.2 is used to implement the interface in the FHD simulations; several comparisons are shown in Section IV.3 to validate the Langevin method of Section III.2.3. The effect of transpiration is demonstrated with several different effusion probabilities, $f$, by examining the long range spatial correlations in the fluid. This is followed by an analysis of the relative computational efficiency of the DSMC and FHD simulations.

We consider a quasi-one dimensional domain with a gradient induced in the gas by setting different temperatures at $x=0$ and $x=L$ (see Fig. 1). These temperatures are set in the DSMC simulations using a fully diffuse boundary condition, and in the FHD simulations with an isothermal no-slip boundary condition. Rather than using a temperature slip boundary condition, the Dirichlet boundary condition for temperature is set to match the mean extrapolated temperatures from the DSMC simulations – this is discussed further below. For the DSMC simulations the $y$ and $z$ boundaries are taken to be periodic; as discussed above, the FHD simulations are performed in 1D. The transpiration interface is at $x=L/2$, the center of the domain. Simulation parameters are given in Table 1.

IV.1 Temperature profiles and effusion condition

In Fig. 2, the temperature profiles are compared between the DSMC and FHD simulations (using the master equation algorithm of Section III.2.2). Comparisons are shown for a range of $f$ values, and in the absence of a membrane. In all cases the temperature boundary conditions for the DSMC simulations are set to 273K and 519K on the left and right hand boundaries, respectively, giving a large gradient of $7604\,\mathrm{K}/\mathrm{cm}$. Due to this large gradient, the gas at the boundary is unable to completely accommodate to the wall temperature, resulting in a significant Knudsen layerSone (2000) and ‘temperature slip’. This feature is correctly captured by the DSMC simulations, but must be accounted for in the FHD simulations with a modified boundary condition. As mentioned above, the FHD boundary conditions have been extrapolated from these new temperatures; see Table 2.

In all cases the DSMC and FHD simulations show good agreement aside from a small discrepancy near the interface, which becomes prominent at $f\gtrsim 0.1$. As discussed in sections I and III, the transpiration models introduced for the FHD simulations are valid only when the effusion condition holds. Similar to the temperature boundary condition discussed above, this is highlighted by the appearance of a Knudsen layer at the surface of the interface, with the relevant Knudsen number being proportional to the effusion probability $f$. For $f\gtrsim 0.1$, this causes the solution near the interface to be poorly approximated by the FHD simulations. This occurs because in this regime there is a substantial chance of a particle crossing the interface in one direction and being reflected back in the other direction after only a small number of collisions, preventing the regions close to the interface from maintaining local equilibrium. This effect is accurately captured by DSMC, which does not rely on a local equilibrium assumption, but not by the FHD simulations.

IV.2 Measurement of hydrodynamic fluctuations

The fluctuations of the hydrodynamic variables are measured by accumulating statistical samples of the conserved variables. For example, in DSMC the mass and momentum densities in the cell centered at $x$ at time $t$ are,

[TABLE]

where $N$ is the number of DSMC particles in the cell, and $\mathbf{v}_{i}$ is the velocity of particle $i$. In FHD these conserved densities, $\rho$ and $\mathbf{p}=\rho\mathbf{u}$, are computed directly (see Eqns. 1-3). Means and second moments are computed from these statistical samples, for example,

[TABLE]

and

[TABLE]

where $S$ is the number of time samples. Correlations of hydrodynamic fluctuations are obtained from these conserved variables, e.g.,

[TABLE]

The calculation of temperature fluctuations is similar; see Ref. 49 for details. For all the results shown in this article the simulation was run until a statistical steady state was reached, recording of the quantities of interest was then begun. The value $S$ is given in the caption for each result. In practice each statistical sample was constructed from an ensemble of shorter simulations to give a total of $S$ time samples.

IV.3 Hydrodynamic fluctuations in DSMC and FHD

In Fig. 3 long range spatial correlations, in the absence of an interface, are compared between DSMC and FHD simulations. Three correlations are shown, $\langle\delta u_{x}^{*}\delta\rho\rangle$, $\langle\delta T^{*}\delta\rho\rangle$, and $\langle\delta T^{*}\delta T\rangle$, in (a), (b), and (c), respectively. In each case the variable with the asterisk is measured at the fixed location $x^{*}=L/4-\Delta x/2$, i.e., in the cell near the center of the left hand side of the domain; see Fig. 1. Note that in the absence of a temperature gradient (i.e., at thermodynamic equilibrium) $\langle\delta u_{x}^{*}\delta\rho\rangle=\langle\delta T^{*}\delta\rho\rangle=0$ and $\langle\delta T^{*}\delta T\rangle=0$ for $x\neq x^{*}$.

We see good agreement between the FHD and DSMC results, observing that the FHD simulations have difficulty capturing the sharp peak in $\langle\delta u_{x}^{*}\delta\rho\rangle$; this confirms the observations of Ref. 19. In (c), the expected Kroenecker delta peak at $x=x^{*}$ has been removed, from both the FHD and DSMC results, for ease of viewing. The FHD simulations again have difficulty in this case, and exhibit significant undershoot in the two adjacent cells. This is a well known effect,Garcia et al. (1987) and these data points have also been removed from the FHD results for ease of viewing – they are shown on the inset plot.

In Figs. 4 and 5, the same correlations are shown in the presence of interfaces with effusion probability $f=0.1$ and $f=0.05$, respectively. In the FHD simulations, the master equation method of Section III.2.2 has been used to simulate the interface. For the $f=0.05$ case, subject to the same caveats as above, good agreement is observed between the FHD and DSMC results in most of the domain, with some scatter in the FHD results in the vicinity of the interface. In the $f=0.1$ comparison, a small systematic difference is observed between the two sets of results, particularly visible in Figs. 4(a) and 4(c) – this represents the breakdown of the effusion condition.

In Fig. 6 the master equation and Langevin approaches for modeling the interface, given in sections III.2.2 and III.2.3, are compared. The correlations $\langle\delta T^{*}\delta\rho\rangle$ and $\langle\delta T^{*}\delta T\rangle$ are shown for the cases $f=0.1$ and $f=0.00625$. Excellent agreement is observed, demonstrating the equivalence of the two approaches.

As can be seen from Figs. 3-5, the general impact of the interface is to reduce the magnitude of long range correlations in the gas. The value $\langle\delta u_{x}^{*}\delta\rho\rangle$ approaches zero for $x/L>0.5$ because the interface does not conserve momentum. Interestingly, a new peak in $\langle\delta T^{*}\delta T\rangle$ develops immediately to the left of the membrane – this behavior is borne out for lower values of $f$ not shown here. This distortion is also observed in a simple Fourier model, with the interface modeled as a region of low thermal conductivity.

Aside from these features, the reduction in magnitude appears to arise principally from the change in temperature gradient. In Fig. 2, a reduction in the gradient is visible with decreasing $f$, where the interface behaves increasingly like an adiabatic surface. In Fig. 7, the magnitude of correlations are shown as a function of the gradient. The top row shows measurements of the value $\overline{|\langle\delta T^{*}\delta\rho\rangle|}$, where the overline indicates the value has been smoothed by averaging over a spatial interval (see figure caption). Fig. 7(a) shows measurements from the left hand side of the interface, while Fig. 7(b) shows the values from the right. In the bottom row, Figs. 7(c) and Fig. 7(d) show the left and right side measurements of $\overline{|\langle\delta T^{*}\delta T\rangle|}$. The results displayed are taken from DSMC and FHD simulations in the range $0.00625\leq f\leq 0.1$, and the subsets shown in each plot are those where a statistically converged solution was obtained. Where both DSMC and FHD results are shown, the DSMC solution should be taken to be the more accurate; we have included the FHD results to illustrate the divergence of the solutions due to the violation of the effusion condition. For lower values of $f$, the magnitude of the correlations becomes smaller relative to equilibrium noise, requiring a larger sample size to obtain a result. FHD results were obtainable in this regime because this approach exhibits greater computational speed; this is discussed further in Section IV.4.

In Figs. 7(a) and 7(b), $\langle\delta T^{*}\delta\rho\rangle$ exhibits a linear response with respect to gradient; this agrees with the findings of Refs. 41 and 51, where a linear response to the temperature gradient is observed in the absence of a membrane. In the same references, a quadratic response is shown in $\langle\delta T^{*}\delta T\rangle$. This appears to be supported by Fig. 7(d), in Fig. 7(c) measurement of this effect is confounded by the appearance of the additional peak near the surface of the interface, discussed above. Finally, it is clear from Figs. 7(b) and 7(d) that the fluctuations on opposite sides of the interface remain correlated even for small values of $f$.

IV.4 Comparative performance

The cell size used to make the comparisons shown in Figs. 3-5, given in Table 1, has been selected based on the kinetic length scale that constrains the DSMC simulations. The time step has been selected to minimize the chance of fluctuations causing negative density or energy values in the FHD simulations - these may also be avoided by increasing the cell size, which is not bound by the kinetic scale as in the DSMC simulations.

In Fig. 8, the effect of changing cell size and time step is compared between the two methods, using the case $f=0.1$. To maintain a constant diffusive stability value, $C\sim\Delta t/\Delta x^{2}$, for every doubling of the cell size the time step is increased by a factor of four. Due to the very small time steps used, the error appears to be dominated by the cell size, with both methods showing the expected quadratic convergence. The FHD algorithm is deterministically third order in time and second order in space and weakly second order for stochastic systems; DSMC is at best second order. We note also that, depending on the simulation parameters, one FHD time step executes approximately ten to thirty times faster than the equivalent DSMC time step.

V Concluding Remarks

In this article, we have investigated the effect of a simple transpiration interface on long range hydrodynamic correlations in a dilute gas, where a non-equilibrium steady state has been induced by a temperature gradient. Several qualitative differences in the correlations were observed, however the principal effect of the interface was to reduce their magnitude proportionally to the effusion probability, $f$. This outcome can largely be explained by changes in the temperature gradient induced by the interface. We also observe a distortion of the $\langle\delta T^{*}\delta T\rangle$ correlations, specifically the appearance of a new peak located near the interface.

The investigation used two distinct numerical techniques, DSMC and FHD, which were found to be in agreement for appropriate values of $f$; for $f\gtrsim 0.1$, non-continuum effects reduce the accuracy of the FHD method. The FHD simulations used both master equation and Langevin approaches for modeling the interface, with the two methods producing equivalent results. This validates the use of FHD for interfaces with low effusion probabilities, giving a computationally faster alternative to DSMC.

There are a number of potential extensions to the work presented here. Currently, both methods used to model the interface in the FHD simulations assume $J_{\mathbfcal{P}}=0$, i.e., there is zero momentum flux. Future extensions to the interface model could include momentum fluctuations around a zero mean, a more realistic model of thermal processes and incorporation of multicomponent effects such as species dependent permeability.

These developments would provide the starting point for numerical modeling of a range of membrane behaviors. More complex models could include models for transport of ions through membranes and the incorporation of chemical reactions that enable active transport. These capabilities would enable the methodology be used to model transport in biological membranes and transport process important to artificial photosynthesis.

Acknowledgements.

The authors wish to acknowledge many years of fruitful discussions with Graeme Bird. This work was supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Applied Mathematics Program under contract DE-AC02-05CH11231. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.