On the Suppression and Distortion of Non-Equilibrium Fluctuations by Transpiration
Daniel R. Ladiges Andrew J. Nonaka John B. Bell Alejandro L. Garcia







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















Center for Computational Sciences and Engineering,
Lawrence Berkeley National Laboratory






Dept. Physics and Astronomy, San Jose State University


 


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.





[
    Alejandro L. Garcia
    December 30, 2023
=======================






§ 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 engineeering<cit.> to molecular biology<cit.>. 

Fluctuations also yield information about the transport properties of the fluid itself<cit.>, and are important in dynamic processes such as phase transitions, combustion and ignition, nucleation, hydrodynamic instabilities, and a host of Brownian dynamics phenomena.<cit.>


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 dynamics<cit.> 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 Carlo<cit.> (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 equation<cit.> so it can be used for flows of arbitrary Knudsen number,<cit.> 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 <cit.>, 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 scattering<cit.> and shadowgraph experiments<cit.>. Numerical methods for solving the compressible and incompressible FHD equations are discussed in Refs. , and the approach has been further extended to multispecies mixtures,<cit.> multiphase flows,<cit.> reacting flows,<cit.> and electrokinetic flows.<cit.>. Hybrid methods combining FHD and DSMC are investigated in Refs. .

When a fluid is held in a non-equilibrium state, long range correlations in the thermal fluctuations<cit.> can result in macroscopic effects.<cit.>


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 dynamics<cit.> 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 <ref> 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 <ref> 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 <ref>, 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 <ref>.





§ 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 <ref>.



 §.§ 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.<cit.> 
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. , , and  with a complete description given in Ref. .

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.<cit.>
When each DSMC particle represents one molecule in the gas, the magnitude of the hydrodynamic fluctuations is in agreement with equilibrium statistical mechanics.<cit.>
For both equilibrium and non-equilibrium systems DSMC yields the physical spectra of spontaneous fluctuations as predicted by theory<cit.> and observed in experiments<cit.>.
In this paper the numerical measurements from DSMC simulations are used to validate the fluctuating hydrodynamic simulation results.



 §.§ Fluctuating Hydrodynamics (FHD)


The equations of fluctuating hydrodynamics are obtained by adding white noise to each dissipative flux in the deterministic Navier-Stokes equations,<cit.> to represent thermal fluctuations. We consider the compressible fluctuating hydrodynamic equations for a pure monatomic ideal gas; generalizations are found in Refs. .


The continuity, momentum and energy conservation equations for FHD are


    ∂/∂ t( ρ)  = - ∇·( ρ u),



    ∂/∂ t( ρ u) =  - ∇·( ρ u ⊗ u ) - ∇ P - ∇·[ Π + Π],



    ∂/∂ t( ρ E )  = - ∇·( ρ u E + p  u) - ∇·[  Q + Q] -
    ∇·( [ Π + Π] ·  u).


Here, ρ is the mass density of the gas, u is the velocity, P is the pressure, T is temperature, and E = c_V T + 1/2| u|^2 is the total specific energy. For a monatomic ideal gas the specific heat at constant volume is c_V = 3/2 k_B /m, where k_B is Boltzmann's constant and m is the molecular mass; the equation of state is P = n k_B T where n=ρ/m is the number density.
From Eqs. (<ref>)-(<ref>) we define the fluxes in the mass, momentum, and total energy as J_M, J_P, and J_ℰ, respectively.

In the momentum and energy equations (Eqs. (<ref>) and  (<ref>)), the deterministic viscous stress tensor Π and deterministic heat flux vector Q are given by 


    Π_ij = -η( ∂ u_i/∂ x_j + ∂ u_j/∂ x_i) + δ_ij( 2/3η∇· u),
        and     Q =  -  κ∇ T,

where δ_ij is the Kronecker delta, η(T) is the shear viscosity, and κ(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. (<ref>) and (<ref>) denote stochastic contributions from Gaussian random fields. The stochastic viscous tensor, Π, stochastic heat flux tensor, 𝐐, 
have zero mean and covariances <cit.>


    ⟨Π_ij(𝐫,t)   , Π_kl(𝐫',t') ⟩
    = 2 k_B η T  [ ( δ_ikδ_jl+ δ_ilδ_jk )
    - 1/3δ_ijδ_kl ]
     δ(t-t') δ( 𝐫 - 𝐫' ),

and


    ⟨Q_i(𝐫,t)   , Q_j(𝐫',t') ⟩ = 2 k_B κ T^2  δ_ijδ(t-t') δ( 𝐫 - 𝐫' ).


where δ(..) is the Dirac delta function.

The FHD Eqs. (<ref>)-(<ref>) 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. .




§ TRANSPIRATION INTERFACE MODELS


As discussed in Section <ref>, 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 condition<cit.>. 



 §.§ 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, ϵ, is conserved but not its momentum. Therefore the velocity of a transmitted particle, 𝐯, is reset by randomly selecting it from a hemisphere of radius |𝐯| = √(2 ϵ/m). 



 §.§ 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.



  §.§.§ 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 by<cit.>

    τ_e = m/ρ f A√(2π m/k_B T).

By comparison, the mean time between collisions for a hard sphere particle<cit.> is

    τ_c = m/4 √(π)ρ d^2√(2π m/k_B T),

where d is the collisional diameter of the particle. The effusion condition described above is n A λ_cτ_e≫τ_c where λ_c is the mean free path, which is satisfied when f ≪ 1. 
The probability distribution for the kinetic energy, ϵ, of particles passing through the interface is,

    P_ϵ (ϵ) = ϵ/k_B^2 T^2 e^  -ϵ/(k_B T);

this is a gamma distribution, Γ(2,k_B T).<cit.> Eqs. (<ref>) and (<ref>) form the basis of the stochastic simulation algorithm given in Section <ref>.

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

    W_→(ϵ) =  f A   n_L /√(2π m k_B T_L)ϵ/k_B T_L e^  -ϵ/(k_B T_L),


    W_←(ϵ) =  f A   n_R /√(2π m k_B T_R)ϵ/k_B T_R e^  -ϵ/(k_B T_R),

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,ℰ), i.e., the probability density of the state on the right of the interface with mass M and total energy ℰ,

    d/d tP(M, ℰ)=∫_0^∞ P(M-m,ℰ-ϵ)    W_→(ϵ) +  P(M+m,ℰ+ϵ) W_←(ϵ) dϵ
    
    -   ∫_0^∞ P(M,ℰ) W_→(ϵ) +  P(M,ℰ) W_←(ϵ) dϵ.
















Following from the above, <cit.> the average one-way flux of mass passing from left to right through the interface is,

    ⟨ J_M⟩_→ = ρ_L f A √(k_B T_L/2π m),

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

    ⟨ J_M⟩= ⟨ J_M⟩_→ -⟨ J_M⟩_←
    =   f A  k_B^1/2/√(2 π m  )(ρ_L T^1/2_L- ρ_R T^1/2_R  ).

 The average energy of a particle crossing the interface is 2 k_B T, so the average total energy flux is

    ⟨ J_ℰ⟩= 2 k_B /m ( T_L ⟨ J_M⟩_→- T_R ⟨ J_M⟩_←)= 2 f A  k_B^3/2/ m √(2 π m  )(ρ_L T^3/2_L - ρ_R T^3/2_R  ).

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 <ref>).

The variances and covariance of the mass and energy fluxes are

    ⟨ J_δ Mδ M⟩ =  m f A  k_B^1/2/√(2 π m  )(ρ_R T^1/2_R+ρ_L T^1/2_L  ),


    ⟨ J_δℰδℰ⟩ =  6 f A  k_B^5/2/ m√(2 π m  )(ρ_R T^5/2_R+ρ_L T^5/2_L  ),


    ⟨ J_δℰδ M⟩ =  2 f A  k_B^3/2/√(2 π m )(ρ_R T^3/2_R+ρ_L T^3/2_L  ).

The resulting Langevin equations are given by

    d /dtM_R= ⟨ J_M⟩ + J_M,


    d /dtℰ_R= ⟨ J_ℰ⟩+J_ℰ,


    d /dtM_L=  -d /dt M_R,
      d /dtℰ_L=  -d /dtℰ_R,

where J_M and J_ℰ are Gaussian noise terms, with

    ⟨J_M(t)J_M(t')⟩=⟨ J_δ Mδ M⟩δ(t-t'),

    ⟨J_ℰ(t)J_ℰ(t')⟩=⟨ J_δℰδℰ⟩δ(t-t'),

    ⟨J_ℰ(t)J_M(t')⟩=⟨ J_δℰδ M⟩δ(t-t').

These equations form the basis of the Langevin method discussed in Section <ref>.






  §.§.§ Master equation implementation

Following from the theory presented above, we define a stochastic simulation algorithm (SSA)<cit.> 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, Δ 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 ρ_L and ρ_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. (<ref>) the rate of particles crossing from each direction is r_L = 1/τ_e(ρ_L,T_L) and r_R = 1/τ_e(ρ_R,T_R).

  * Initialize the accumulation of mass and energy flux, F_M and F_ℰ, across the interface to zero.

  * Set the interface time counter t_ssa = 0. Loop while t_ssa < Δ t:


  * Given an exponential distribution of time

between particle crossings, the time until the next particle crossing is given by

    θ = -log(ℛ_1)/ r_L + r_R,

where ℛ_i ∈ (0,1] indicates an independent uniform random number.

  * Increment the interface time counter, t_ssa+= θ. If t_ssa≥Δ t then break out of the loop.

  * Select the direction for the particle crossing the interface. If r_L/(r_L + r_R) > ℛ_2, then the particle crosses from left to right, otherwise it crosses from right to left.

  * Following from Eq. (<ref>), the energy of the effusing particle is gamma distributed and generated by<cit.>

    ϵ = -k_B T_* log(ℛ_3ℛ_4),

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 += m and F_ℰ+= ϵ. Otherwise F_M -= m and F_ℰ-= ϵ.


  * Finally, the left cell is updated as ρ_L -= F_M /Δ V and ρ_LE_L -= F_ℰ /Δ V, and the right cell is updated as ρ_R += F_M /Δ V and ρ_R E_R += F_ℰ /Δ V with Δ 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.



  §.§.§ 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 ρ_L and ρ_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. (<ref>) and (<ref>), and the variances and covariance given by  Eqs. (<ref>)-(<ref>).

  * Calculate the correlation coefficient, given by

    r=⟨ J_δℰδ M⟩/√(⟨  J_δℰδℰ⟩⟨  J_δ Mδ M⟩).


  * Generate correlated random variables,

    𝒢_3 = 𝒢_1 r + √(1-r^2)𝒢_2

where 𝒢_1 and 𝒢_2 are independent Gaussian random numbers of mean zero and variance one.

  * Following from Eq. (<ref>), the mass and energy flux accumulations are given by

    F_M = Δ t ⟨  J_M⟩ + √(Δ t⟨  J_δ Mδ M⟩) 𝒢_1,


    F_ℰ = Δ t ⟨  J_ℰ⟩ + √(Δ t⟨  J_δℰδℰ⟩) 𝒢_3.


  * Finally, the left cell is updated as ρ_L -= F_M /Δ V and ρ_LE_L -= F_ℰ /Δ V, and the right cell is updated as ρ_R += F_M /Δ V and ρ_R E_R += F_ℰ /Δ V with Δ 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.




§ 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 <ref> is used to implement the interface in the FHD simulations; several comparisons are shown in Section <ref> to validate the Langevin method of Section <ref>. 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. <ref>). 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 <ref>.










 §.§ Temperature profiles and effusion condition



In Fig. <ref>, the temperature profiles are compared between the DSMC and FHD simulations (using the master equation algorithm of Section <ref>). 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 K/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 layer<cit.> 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 <ref>. 

In all cases the DSMC and FHD simulations show good agreement aside from a small discrepancy near the interface, which becomes prominent at f ≳ 0.1. As discussed in sections <ref> and <ref>, 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 ≳ 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. 











 §.§ 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,

    ρ(x,t) = m/Δ V N(x,t)
        and    𝐩(x,t) = m/Δ V∑_i^N 𝐯_i,

where N is the number of DSMC particles in the cell, and 𝐯_i is the velocity of particle i. In FHD these conserved densities, ρ and 𝐩 = ρ𝐮, 
are computed directly (see Eqns. <ref>-<ref>). Means and second moments are computed from these statistical samples, for example,

    ⟨𝐩(x^*) ρ(x) ⟩ 
    = 1/S∑_t^S 𝐩(x^*,t) ρ(x,t)

and

    ⟨δ𝐩^* δρ⟩ =
    ⟨δ𝐩(x^*) δρ(x) ⟩ =
    ⟨𝐩(x^*) ρ(x) ⟩
    - ⟨𝐩(x^*) ⟩⟨ρ(x) ⟩,

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

    ⟨δ𝐮^* δρ⟩ = 
    1/⟨ρ^* ⟩( 
    ⟨δ𝐩^* δρ⟩
    - ⟨𝐮^* ⟩⟨δρ^* δρ⟩).

The calculation of temperature fluctuations is similar; see Ref.  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.



 §.§ Hydrodynamic fluctuations in DSMC and FHD



In Fig. <ref> long range spatial correlations, in the absence of an interface, are compared between DSMC and FHD simulations. Three correlations are shown, ⟨δ u_x^* δρ⟩, ⟨δ T^* δρ⟩, and ⟨δ T^* δ T ⟩, in (a), (b), and (c), respectively. In each case the variable with the asterisk is measured at the fixed location x^* = L/4-Δ x/2, i.e., in the cell near the center of the left hand side of the domain; see Fig. <ref>. Note that in the absence of a temperature gradient (i.e., at thermodynamic equilibrium) ⟨δ u_x^* δρ⟩ = ⟨δ T^* δρ⟩ = 0 and ⟨δ T^* δ T ⟩ = 0 for x ≠ x^*.

We see good agreement between the FHD and DSMC results, observing that the FHD simulations have difficulty capturing the sharp peak in ⟨δ u_x^* δρ⟩; this confirms the observations of Ref. . 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,<cit.> 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. <ref> and <ref>, 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 <ref> 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. <ref> and <ref> – this represents the breakdown of the effusion condition.







 

In Fig. <ref> the master equation and Langevin approaches for modeling the interface, given in sections <ref> and <ref>, are compared. The correlations ⟨δ T^* δρ⟩ and ⟨δ T^* δ T ⟩ 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. <ref>-<ref>, the general impact of the interface is to reduce the magnitude of long range correlations in the gas. The value ⟨δ u_x^* δρ⟩ approaches zero for x/L > 0.5 because the interface does not conserve momentum. Interestingly, a new peak in ⟨δ T^* δ T ⟩ 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. <ref>, a reduction in the gradient is visible with decreasing f, where the interface behaves increasingly like an adiabatic surface. In Fig. <ref>, the magnitude of correlations are shown as a function of the gradient. The top row shows measurements of the value |⟨δ T^* δρ⟩|, where the overline indicates the value has been smoothed by averaging over a spatial interval (see figure caption).
Fig. <ref> shows measurements from the left hand side of the interface, while Fig. <ref> shows the values from the right. In the bottom row, Figs. <ref> and Fig. <ref> show the left and right side measurements of |⟨δ T^* δ T ⟩|. The results displayed are taken from DSMC and FHD simulations in the range 0.00625 ≤ f ≤ 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 <ref>.

In Figs. <ref> and <ref>, ⟨δ T^* δρ⟩ exhibits a linear response with respect to gradient; this agrees with the findings of Refs.  and , 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 ⟨δ T^* δ T ⟩. This appears to be supported by Fig. <ref>, in Fig. <ref> 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. <ref> and <ref> that the fluctuations on opposite sides of the interface remain correlated even for small values of f.






























 §.§ Comparative performance

The cell size used to make the comparisons shown in Figs. <ref>-<ref>, given in Table <ref>, 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. <ref>, 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∼Δ t / Δ 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.







§ 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 ⟨δ T^* δ T ⟩ 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 ≳ 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_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.



































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.





