Transmission probability of gas molecules through porous layers at Knudsen diffusion

Gas flow through layers of porous materials plays a crucial role in technical applications, geology, petrochemistry, and space sciences (e.g., fuel cells, catalysis, shale gas production, and outgassing of volatiles from comets). In many applications the Knudsen regime is predominant, where the pore size is small compared to the mean free path between intermolecular collisions. In this context common parameters to describe the gas percolation through layers of porous media are the probability of gas molecule transmission and the Knudsen diffusion coefficient of the medium. We show how probabilistic considerations on layer partitions lead to the analytical description of the permeability of a porous medium to gas flow as a function of layer thickness. The derivations are made on the preconditions that the molecule reflection at pore surfaces is diffuse and that the pore structure is homogenous on a scale much larger than the pore size. By applying a bi-hemispherical Maxwell distribution, relations between the layer transmission probability, the half-transmission thickness, and the Knudsen diffusion coefficient are obtained. For packings of spheres, expressions of these parameters in terms of porosity and grain size are derived and compared with former standard models. A verification of the derived equations is given by means of numerical simulations, also providing evidence that our analytical model for sphere packing is more accurate than the former classical models.


Introduction
Models of gas flow through porous media can be categorized according to the assumptions made on the properties of the solid material and the gas species flowing through the pores.For inert gases the models essentially differ in the way the pore structure and the gas collisions are represented.A systematic approach to the description of rarefied gas flow through porous media dates back to Knudsen's work on gas flow through cylinders [1].The often used application to porous media regard the pores in the medium as a collection of straight cylindrical filaments which are distributed arbitrarily (with varying orientations) throughout the medium.More sophisticated models were considered later by Derjaguin [2,3] and Asaeda et al. [4], which are based on packing of spheres.These approaches are also used in the present article to link the transmission probability with other parameters used to describe the gas diffusion through the porous medium.Recent works apply numerical methods, in particular Direct Simulation Monte Carlo (DSMC) [5,6], to take complicated pore shapes into account [7,8].
Here we focus on the gas density range where the pore size is much smaller than the mean free path between intermolecular collisions.This Knudsen regime plays an important role in material sciences, especially in relation with adsorption [9,10] and catalytic processes [11], which are the basis of many technical applications (e.g., fuel cells and catalysts for a great variety of chemical products).It was studied also in the context of shale gas extraction [12], which became more important for the worldwide energy supply in the last decades.Furthermore, the Knudsen pressure regime is of great relevance in space sciences.In particular, it prevails at processes of outgassing from comet surfaces [13][14][15].At sufficiently low pressures the so-called Knudsen diffusion determines the gas flow.In contrast to viscous flow, where the pressure gradient causes the upstream molecules to impart, on the average, a momentum to the downstream molecules, the situation is different for Knudsen diffusion: Because the molecule motion in the pore volume occurs without intermolecular collisions, a net flow can only be the result of spatial variations of gas density or temperature (mean molecule speed).Assuming diffuse or specular reflection of gas molecules at the pore walls, the Knudsen diffusion through a layer of porous material depends only on the pore structure and the density/velocity distribution of the gas molecules.Since the motion of gas molecules through the medium is a stochastic process, the steady-state flow through the layer can be analyzed on the basis of probabilistic considerations, as elaborated in this article for the special case of planar layers.Among other relations it yields the fraction T (z) of molecules reaching the distance z from the surface in the porous medium (the others return to the surface before having reached this depth).T (z) can also be regarded as the molecule transmission probability of a layer of thickness z, i.e., the probability that a molecule entering the pores at one side of the layer will exit at the other side.Former authors like Gundlach et al. [13], Skorov and Blum [16] used the heuristic formula with constant medium parameters a and b.The probability approach described in the following enables us a relatively simple justification of this formula.For this purpose the scenario illustrated in Fig. 1 is investigated, which consists of an infinite planar layer of porous medium divided into two sub-layers.Most of the molecules entering the layer at one side (boundary B 0 ) leave the medium only after long erratic paths, since they are reflected many times at the pore wall in the medium.The fraction of molecules penetrating to a certain distance z from the surface B 0 decreases with increasing z.Therefore, when the inflow at the opposite sides of the layer is different, a corresponding change of the gas density prevails across the layer under steady-state conditions.Neglecting intermolecular collisions and providing isothermal conditions (implying constant mean molecular speed across the sample), the net molecular flux divided by the magnitude of the gas density gradient is a constant D K , the Knudsen diffusion coefficient.
In many applications it is advantageous to use T to describe the gas molecules throughput through a layer [13,14].In other circumstances, in particular when the gas flux or density (as it develops throughout the medium) is of interest, the Knudsen diffusion coefficient D K is sought [15,17,18].One of the possible approaches to link these parameters is based on a bi-hemispheric Maxwell model of the velocity distribution, which was also applied by Asaeda et al. [4] for gas flow through packings of spherical grains.They expressed D K as a function of the grain diameter d g and the porosity of the packing.Further details are discussed in Sect. 4.
Section 2 starts with an investigation of the transmission properties of two layers of constant porosity, which is related to the Knudsen diffusion constant in Sect.3, and compared with models for D K in Sect. 4. Section 5 gives a validation of the derived analytical models by means of numerical simulations on the basis of two independent approaches.The most important symbols used in this article are illustrated in the Figs. 1 and 2 or listed in Table A1 for convenience.

Transmission through two layers
The analysis of the interdependence of transmission probability and diffusion coefficient requires some basic relations on the transmission and reflection properties of the porous layer.For this purpose the medium is divided into two sub-layers of thickness L 1 and L 2 as shown in Fig. 1, with outer boundaries B 0 and B 2 and interfacial boundary B 1 .We assume that the medium in the layers is homogeneous, with porosity .Each boundary is considered as a planar medium cross section.For a packing of spherical grains of diameters d g it means that the spherical sections projecting beyond B 0 and B 2 are removed.This plane surface is different from surfaces of real-world packings where the surface consists of spherical surface sections.However, the plane-cut assumption facilitates a clear distinction between the immediate reflection of incident gas molecules at the surface boundary and the reflection by their coming back to the surface after having entered the layer(s) pores.We further assume that the cross-section area porosity (area of pore cross section divided by total cross-section area) equals the volume porosity of the whole layer.This precondition is guaranteed by the main assumption maintained throughout this article, which is the homogeneous distribution of grains in each sub-layer.More precisely, the distribution of the orientations of pore surface area elements does not depend on location (x, y, z) on a macroscopic scale much greater than the typical grain/pore size.
As a gas molecule traverses two adjacent layers or is reflected by them, it may cross the boundary between the layers several times before exiting the medium on either side.Possible path examples are indicated in Fig. 1: (a) Immediate reflection at the surface B 0 , which occurs with probability 1 − , provided the above-mentioned preconditions.(b) The gas molecule enters the layer 1, moves through its pores, and comes out after one or several reflections at the pore walls.(c) The molecule also enters layer 2 but escapes again through the surface B 0 after having moved forth and back between the layers, crossing the boundary B 1 an even number of times.(d)&(e) The molecule does finally leave through boundary B 2 , which may happen after crossing B 1 one or several (an odd number of) times.The cases (a), (b), and (c) are reflections at the two-layer structure, whereas (d) and (e) are transmission scenarios which may be much more involved than shown in Fig. 1 (many crossings of B 1 are possible).Let T i be the transmission probability of layer i, and R i the corresponding layer reflection probability, under the assumption that the molecule has entered the layer pores (has not been reflected at the grain cut areas of the entrance surface B 0 ).We consider an inert solid matrix without adsorption processes at the pore surface.So T i + R i = 1 since no sinks or sources are present.The probability of a molecule entering the pores through the boundary B 0 equals .So the probabilities of the paths (a)-(e) are 1 respectively.The probability of a molecule incident onto B 0 to traverse both layers after moving forth and back between the layers, crossing B 1 (2n + 1) times, is T 1 (R 2 R 1 ) n T 2 .The total transmission probability T t is the sum over these possibilities Here T denotes the transmission probability under the condition that the molecule crosses the boundary at a pore opening.The probability that the molecule is reflected by the two-layer structure can be calculated by a similar sum of path probabilities, or simply as where R is the reflection probability associated with the two-sub-layer structure under the condition that the molecule falls onto an open pore.Dividing reflection by transmission probability proves that this fraction is an additive quantity Thus, by dividing the layer in more and more sub-layers (with equal properties), we recognize that R/T is proportional to the layer thickness L. So it can be interpreted as a measure of resistance of the layer to gas flow, and R/(T L) must be constant (independent of L).In other words, the layer reflection and transmission probabilities, R(L) and T (L), which are functions of the layer thickness L, can be used to define the quantity which is independent of L and characterizes the medium itself.By substituting T = 1 − R in (7) we obtain after rearrangement These equations reveal that L h is the layer thickness for which half the gas molecules incident on open pores traverse the layer through the pores, the other half returning to the surface boundary where they entered.The last result can be rephrased: among the gas molecules entering the layer, the fraction 1/(1 + z/L h ) reaches the distance z from the surface.It does not decrease exponentially with z, which is a result of the fact that molecules reflected backwards are not removed from the flow process (as it occurs for molecule and light ray attenuation due to absorption), but may come again to the same position z after other reflections at the pore wall and then may advance further into the medium.The total transmission probability T t is obtained from ( 9) by taking into account (2).It is, in fact, of the form (1) discussed above, showing the meaning of a = and b = L h .Since in reality the boundary of a sample is not a plane cut, the parameters a and b may differ from these values, allowing for a modification of the reflection at the layer boundary in accordance with its uneven surface structure.

Relation between transmission property and diffusion coefficient
A relation between L h and D K can be found by explicitly analyzing the flow contributions of the gas molecules moving forwards and backwards (positive and negative velocity components v z , respectively), giving a net forward flux, i.e., in +z direction, normal to the layer surfaces.In Fig. 1 this forward direction is vertically upwards, with z being the height above B 0 .We assume steady-state isothermal conditions, gas and porous medium having the same temperature T .Provided that the gas is in thermal equilibrium above and below the layer at sufficient distance, the molecule velocities satisfy a Maxwell distribution with average molecule speed where R = 8.31 J/mol K is the ideal gas constant and M the molar mass of the gas species.Near the lower surface, the molecules moving towards the layers satisfy the corresponding hemispherical Maxwell distribution (v z > 0), but the leaving molecules may have a different distribution.In particular they must have a different density if a steady-state net gas flow is to be maintained.However, it is plausible to assume that also the molecules moving backwards (v z < 0) have an overall hemispherical Maxwell distribution, albeit with a different density.The reason is that the preconditions of an isothermal process and diffuse reflection at the pore surface, with the surface area elements possessing an isotropic distribution of orientations, have these implications: (i) the energy distribution of the reflected molecules is the same as for the incident molecules, resulting in a Maxwell distribution of the molecule speed with the same average c in every direction; (ii) since the molecules 'forget' their direction of incidence at diffuse surface reflection, no direction is predominant among the reflected molecules when the forwards and backwards moving densities are not too different.In this case the vectorial velocity distribution of the molecules returning from the layer is practically Maxwellian.A similar argument holds for molecules entering the other side of the layer.Summing up, we apply a model of two hemispherical Maxwell velocity distributions at each cross section z, with two number densities n + (z) and n − (z) inside the pores for molecules with v z > 0 (moving forwards) and v z < 0 (moving backwards), respectively.This model was already successfully applied by Asaeda et al. [4] to derive D K for a monodisperse packing of spheres, the consequences of which are further studied in Sect. 4.
The partial densities n + (z) and n − (z) are defined as double the actual densities of the molecules moving in one direction, so that we can apply the familiar formula J ± p = n ± c/4 for the molecular flux J ± p inside the pores in the ±z direction (it amounts to the Hertz-Knudsen equation [19,20] by means of (10) and the ideal gas law).The total number density is the average With the help of the split distribution n ± (z) we can relate the transmission probability T and the half-transmission length L h to the Knudsen diffusion coefficient D K .For symmetry reasons the total net flux J is parallel to the z-axis.It is composed of the unidirectional contributions J ± of molecules moving forwards and backwards (in +z and −z direction, respectively) The factor accounts for the fact that inside the medium the gas can only move through the pores, wherefore the flux is reduced by the porosity.Here J is defined as the molar flow rate per cross section of the whole medium, whereas J p = J / is the actual flux through the cross section of the void (open pores) in the medium.The vectorial flux can be expressed as a product of D K and the density gradient ∇n (Fick's law for Knudsen diffusion), so we find the equality In the following the denotation n ± 0 = n ± (0) and n 0 = n(0) is used for brevity, and similarly for z = L. Thus, n 0 and n L are the total gas densities in the pores at the opposite surfaces of a layer ranging from z = 0 to z = L. Utilizing the fact that J (z) is constant for a steady-state flow (in a vessel we have to assume constant cross-sectional area to ensure this condition), the last equation can be integrated over the layer height L to give Under steady-state conditions, the fluxes of the molecules leaving the medium can be expressed in terms of the fluxes of the incident molecules and the probabilities of transmission T and reflection R = 1 − T in the pores.Due to the proportionality (12) the corresponding relations can be written in terms of the respective densities Substitution of these relations into the density-dependent terms of ( 13) and ( 15) yields n + − n − is independent of z due to the steady state of the flow, and the constant flux J can be eliminated from the above equations to obtain According to the definition (7) we finally arrive at the sought expression for the halftransmission thickness These last two relations enable us to determine the transmission properties from models of the Knudsen diffusion coefficient.The most important examples are given in the following section.

Expressions by models of Knudsen diffusion
Three well-known models for the Knudsen diffusion constant D K are used in this section to express the transmission properties in terms of the layer pore structure.The models were developed by Knudsen [1] (see also Epstein [21] for application to porous media), Derjaguin [3] and Asaeda et al. [4], respectively.We add another model for layers of porous media with homogeneous and isotropic distributions of pore surface elements, by deriving the layer reflection probability and applying (20) to obtain D K .In Sect. 5 the good concurrence with numerical simulations will prove it more accurate for sphere packings than the former models.

Knudsen cylinder model
The first model is based on a representation of the pores as a set of slanted cylindrical tubes, which may have an average inclination with regard to the net flow direction.The model is based on Knudsen's formula for the flow through a single cylinder of diameter d c [1].Some adaptations are needed for porous media.Since the volume of the bundle of cylinders occupies only the void fraction , the flux is reduced by the same factor.In addition, an inclination of the cylinders by an average angle θ results in a reduction of the pressure gradient along the cylinders by a factor 1/τ = cos θ as well as an increase of the cylinder cross sections parallel to the layer planes by the factor τ = 1/ cos θ .Both effects decrease the net flux by the same factor 1/τ .Thus, the diffusion coefficient of a bent cylinder bundle becomes [21] The so-called tortuosity τ is typically between 1 and 3, but may also reach higher values.In a homogeneous and isotropic distribution of pieces of cylinder pores throughout the medium one would expect τ = √ 3, but much larger values may occur in compactly packed beds [11], in particular when the grains are elongated or flat with the widest extensions lying normal to the net flow direction.The cylinder model shows the dependence on the cylindrical pore diameter and mean speed of gas molecules, but it does not fully account for the pore structure and the porosity of the medium (since the tortuosity depends on porosity as well).

Derjaguin/Asaeda models
To overcome some of the disadvantages of the cylinder model Derjaguin [2] considered a packing of monodisperse spheres (republished as Derjaguin [3]), the results of which are shortly summarized in the following for its importance.Tracing the gas molecules gives a series of straight chords (rectilinear segments) between successive collisions with sphere surfaces.Adding chords and taking the average change of directions between successive segments into account, Derjaguin determined the diffusion coefficient in analogy to Brownian motion via the average square of the vectorial displacement d gas molecules attain in time t, Here and in the following the bar (overline) expresses the average of the respective quantity (e.g., λ 2 is the average of the squared chord length, whereas λ 2 is the square of the average chord length).In contrast to Derjaguin [3] we define the gas density as the moles per pore volume V p (and not per medium volume), wherefore the denominator appears on the left side of (23).This formula has been obtained under the assumption of diffuse reflection, for specular reflection the result is similar, with the second (subtracted) term in the parentheses omitted.The chord length distribution plays an essential role.According to (23) the first and second moments of the chord length distribution are sufficient to determine the diffusion coefficient for the Knudsen regime, provided diffuse or specular reflection occurs at the pore walls.A further simplification can be achieved by assuming that the molecule collision probability is independent of its past (Markov process), with dλ/λ being the probability of collision in the infinitesimal path element dλ, giving a homogeneous Poisson process with the mean chord length λ between successive collision events.With this assumption the chord lengths λ obey the exponential distribution density exp(−λ/λ)/λ, with λ 2 = 2λ 2 .Assuming that the distribution of chords is homogeneous throughout the pore volume V p and the collisions are distributed uniformly over the pore surface S p , the introduction of the specific surface s p = S p /V p of the pores is possible by means of [3] with porosity and grain (sphere) diameter d g .Substitution into (23) enables us to express D K in terms of and d g with the factor = 13/6 as introduced by Asaeda et al. [4].Here the specific surface of the pores, s p , has been related to the specific surface of the solid, s sol , via A comparison of ( 25) with the Knudsen model shows the dependence of the tortuosity τ on the porosity of the medium.The expression (25) was derived also by Asaeda et al. [4], but on the basis of a completely different approach.As mentioned in the former section a bi-hemispherical Maxwell distribution was considered for this purpose, with different densities for forwards and backwards moving molecules.The momentum loss due to the diffuse reflection on the spherical surfaces was calculated, resulting in the density gradient accompanied with the given flux.The authors arrived at the same equation (last term in (25)), with complicated integrals for the number , which they calculated to be 2.18.However, our numerical calculation of their integrals shows that = 13/6 within quadrature accuracy, so full agreement with Derjaguin's result is verified.

Model based on reflection probability of sub-layers
For a homogenous packing we assume that in each thin sub-layer the distribution of pore surface elements is isotropic (no preferential surface directions).Therefore we calculate the reflection probability of such a sub-layer by summing up the backscattered gas molecules over all pore surface elements, divided by the total number of incident molecules.This is achieved by first considering the fraction of molecules reflected by an infinitesimal element of the pore surface with outward normal n.For this purpose we have to define the reflection at the pore surface in mathematical terms.The diffuse reflection assumption used throughout this article is described by where c is the mean speed of the gas molecules.The total collision rate of forward moving molecules (i.e., with v z > 0) per unit area of the pore surface can be written The main steps in the evaluation of this integral are given in A. As to be expected the collision rate J i of forwards moving molecules (v z > 0) depends only on the colatitude θ = arccos(n • e z ), i.e., the angle between the normal n of the pore surface element and the z-axis, with J i = 0 when the surface normal points in z-direction.
We consider a uniform distribution of pore surface directions in each thin sublayer.The fraction of molecules reflected by the layer is obtained by integrating the reflected molecules over all area elements n d A present in the layer.Due to their isotropic distribution, we can express them by the specific surface s p of the pores.For this purpose we take into account that the pore volume V p of a thin sub-layer can be expressed by its thickness dz, the layer's cross-section area A, and the porosity , The notations d A(n) and d (n) emphasize that the area and solid angle elements are associated with the surface normal direction n. dA(n) emits only a fraction of the incident molecules in directions v z < 0. Since the reflection is diffuse, the outgoing molecules are isotropically distributed in the pore hemisphere beyond the tangential plane at the respective surface point.This means that the surface emits the fraction n • e r d (e r )/π into the solid angle d (e r ).So the rate of molecules reflected backwards at the surface element d A(n) (i.e., molecules incident with v z > 0 and scattered with v z < 0) is where the same integral appears as in (27), but this time integration is over the set (n) of emission directions e r , which are opposites to the set of e i .Substitution of ( 27), (28), and (29), and integration over all surface elements yield the following reflected molecule flow rate: where the integral in the second line is 2/3.Since the incident flow rate of molecules entering the layer from one side with v z > 0 is A cn + /4, the reflection probability of the infinitesimal layer becomes s p dz/3.On the other hand, it equals dz/L h in virtue of equation ( 8) for infinitesimal thickness L = dz.Thus, we obtain a formula for the half-transmission length By means of the relation (21) we finally yield another model for the Knudsen diffusion coefficient This result has been entirely derived on the basis of probabilistic considerations, which is in contrast to Asaeda et al. [4] who applied momentum conservation to the reflection process.The preconditions of the layer model are diffuse reflection at pore surfaces, constant temperature throughout the sample (and gas volume), and a homogeneous and isotropic distribution of pore surface elements in each thin sub-layer on a lateral scale much greater than the sphere diameters.The probability model ( 33) is similar to the formula (25) derived by Asaeda et al. [4], which is a simplification of the more general formula (23) derived by Derjaguin [2].Our result (33) is by a constant factor 13/16 smaller than (25).Since the numerical simulations confirm the higher accuracy of our result (33), we can conclude that the simplification that implies the Asaeda formula (25) from the more general Derjaguin formula (23) does not hold well for random sphere packing.That is, the chord lengths λ do not obey the distribution density exp(−λ/λ)/λ, and λ 2 = 2λ 2 is an estimate that is not quite accurate (as shown in Sect. 5 leading to an error of about 20 % for ideal sphere packings).Table 1 summarizes the relation between the different models when expressed as half-transmission thickness L h .The last line can be obtained from the previous line by the above-discussed approximation of the path length distribution.The equations in Table 1 reveal that the half-transmission length L h is similar to the mean chord length.It is of the same order of magnitude as the equivalent cylinder diameter if the tortuosity is not too high, which depends on the porosity as revealed by comparing (34) with ( 35) and (37).
When equating d c = d g the tortuosity factor alone is used to represent the influence of the grain arrangement (and so the pore structure): Here q = 16/13 or q = 1 when ( 35) or (37), respectively, is used as a reference.The comparison of the Knudsen cylinder model with other models as reference requires different q values.The additional correction factor q can also be chosen in accordance with results from experiments or simulations.The advantage of expressing τ in terms of ( 38) is that it represents well the main porosity dependence for realistic sphere packings (in agreement with the discussed models), leaving only a minor porosity dependence of the correction factor q (an example is given below).According to formula (35), L h agrees with the sphere diameter for packings of porosity ≈ 2/3.A considerable increase of porosity above this value can make the medium much more permeable to gas flow, with a substantial amount of molecules advancing much deeper than a grain size into the layer (L h > d g ).This situation occurs in comet nucleus material where the low gravity allows very high porosities of 70-75%, which may even reach above 80% at certain places [22].
The model (37) by Asaeda et al. [4] was derived under the precondition that all spheres have the same diameter (monodispersity).In contrast, formula (36) by Derjaguin [2] and also the above result (32) are valid for polydisperse packings as well.The Derjaguin model as well as our derivation for the layer model (when expressed in terms of the specific surface s p ) do not assume anything about the grain size and shape, apart from the precondition of an isotropic and homogeneous distribution of surface elements of the pore walls, and that the emission directions of successive collisions are independent of each other.Of course, the assumption of homogeneity amounts to a good mixing of the spheres of all sizes in a polydisperse packing.For this case (35) can be generalized by adding the contributions of all sorts of spheres to the specific surface.So (26) implies where the sum goes over all sorts of spheres with n m being the number density and d m the diameter of the sort m.
In practice, we have to expect that realistic packings (even for spherical grains) are not perfectly homogeneous and depend on the creation process.This holds even more for angular grains.A variation of the spatial distribution of grain diameters can occur due to sample vibration under gravity, causing segregation (Brazil nut effect).Dependent on the packing process, the pore space may include irregularly distributed large cavities or constrictions, and it may contain anisotropic structures on a scale larger than the average chord length of gas molecule motion.In consequence, certain regions may be easier accessible to the gas flow than others, and anisotropy may result in direction-dependent permeability of the medium.Such effects cannot be represented well by the discussed models.An attempt to account for some of these effects is often made by introducing a suitable tortuosity.For instance, Güttler et al. [ [23] have recently shown by DSMC simulations for packings of steel spheres that the introduction of a porosity-dependent correction factor q = 1.6 − 0.73 in (38) substantially improves the agreement of the model with the simulation results.Of course, a location-and direction-dependent character of a sample cannot be represented by this heuristic approach.

Overview of methods
For the assessment of the validity of the formulas ( 21) and ( 35)-(37) we perform numerical simulations for sphere packings of different porosities but constant grain diameters.Since solvers for fluid flow are not applicable in the Knudsen regime, we have to resort to programs which follow the molecules (or group of molecules building a fictive particle) as they travel through the pore space of the medium.Molecular Dynamics (MD) solvers [24] can be used to model the interaction of gas molecules with each other and with the pore surface in detail.Since we focus on diffuse reflection and large number of molecules meandering through a complex pore geometry, MD would take extremely large computation times.The Lattice Boltzmann Method (LBM) can be used in the transition regime to simulate the behavior of a gas molecule ensemble passing the layer [25].However, it is based on a lattice discretization and therefore not the best choice when the possible trajectories of molecules moving through the layer are to be determined.For these reasons we exclude MD and LBM solvers.
The diffuse collisions at the pore wall of the solid matrix are not deterministic but obey a stochastic rule.So the molecule trajectories are random walks, composed of chords between successive collisions with the pore wall.Therefore, the first approach for validation of the analytical models is a specific Random Walk (RW) algorithm.Actually it is a 'Test Particle Monte Carlo' method which has been successfully applied for years in the context of gas transport in cometary material [14].Details on the theoretical background and algorithm can be found in the pertaining literature [26,27], where also further references are provided.In the following we continue to refer to it as Random Walk (RW) due to its nature and the fact that each single particle is followed along its erratic path through the pores.The advantage of this approach (compared to the other applied method) is that it works with sphere packings, handling the spherical geometry precisely, without the need of representing the sphere surface by an appropriate mesh.The RW program is applied to calculate the penetration depth of molecules into the porous layer, giving directly the transmission probability T (z) as a function of the layer thickness z.
As a second numerical technique we apply the Direct Simulation Monte Carlo (DSMC) program SPARTA1 [28].It is used in two ways: first, to verify the penetration depth of gas molecules, and second, to determine the density distribution of the gas flow through the packing under steady-state conditions.The DSMC is dedicated to the simulation of rarified gas flows, based on an algorithm for surrogate molecules, each of which represents a huge number of real molecules.The theoretic justification of the approach and its implementation have been studied for decades.An in-depth description is given by the main developer of DSMC [6].The DSMC can also be applied in the transition regime, where, in addition to Knudsen diffusion, viscous flow contributions start to play a role.However, here we focus on the Knudsen diffusion by switching off intermolecular collisions and only allowing for collisions at the pore surface.In contrast to the applied RW program, the spheres are represented as triangle meshes, which requires more preprocessing of the models.

Packing geometry and preparation
Simulation samples of fixed porosities are generated in a cuboid box, as illustrated in Fig. 3 for the sample with = 0.65.In the simulations the gas molecules enter the box at the bottom (z = 0) and transmitted molecules leave at the top (z = 30 mm) of the box.At the vertical sides periodic boundary conditions are applied.The spheres projecting beyond the sample box are cut off, leaving cuts shown as circular red disks at the box boundary.At the side facing the observer the cuts are omitted for better visibility of the packing.
Special care is taken to prepare packings with random but uniformly distributed spheres.This means that the porosity undergoes only slight changes on a scale much larger than a single sphere.Figure 4 shows the area porosity (the area of the void part divided by the total area of a layer cross section) as a function of the height z in the simulated samples.To avoid significant variation of porosity near the wall of the sample box, the packing is constructed with periodic conditions at opposite box boundaries, so that an infinite layer is well simulated when the simulation domain is about 20 or more sphere diameters wide.All conditions are achieved by the packing routines developed by Vasili Baranov, which he provides via a github repository. 2hree packings of 1 mm spheres are used, with filling factors 0.55, 0.35, and 0.15, covering a wide porosity range ( = 0.45, 0.65, and 0.85).Since the sphere diameter of 1 mm and the size of the simulated sample (cube-shaped box of 30 mm side length) are always the same, different number of spheres need be used to obtain the sought porosities, which amounts to about 30000, 18000, and 8000 for the respective packing.
The following considerations led to the choice of the millimeter sphere diameter: At a fixed gas temperature the studied Knudsen diffusion models can be well transferred to a scaled geometry (where all grain radii and their position vectors are multiplied by the same scale factor, which keeps the porosity constant): The Knudsen diffusion constant D K and the half-transmission thickness L h have to be multiplied by the scale factor, whereas the transmission and reflection probabilities, T and R, stay the same.Therefore the choice of the sphere size for the simulations is not essential when only the molecular diffusion (Knudsen diffusion) is investigated.The results are Fig. 3 Grain packing of 1 mm diameter spheres in a cube of 30 mm side length, as used for the simulation of a layer with porosity 0.65.Sphere cuts at the box boundaries are colored red (for better visibility cuts are omitted at the side facing the observer) Fig. 4 Areal porosity (solid lines) for cross sections at height z in the simulated samples with nominal porosity 0.45, 0.65, and 0.85 (dash lines) applicable, as well, to mesoporous materials or larger grains.The only precondition is that the gas densities stay in the Knudsen regime.For instance, dust on comets covers the whole mentioned size range.Because of the low gas pressure on the comet, the Knudsen regime prevails in pores up to about millimeter size.This is of particular interest in the context of dust agglomerates in near-surface layers, which are of great importance for the thermal balance and outgassing of the comets [30,31].Packing of millimeter size are also very suitable for laboratory measurements to study the Knudsen diffusion because of handling and safety aspects.This holds especially for spherical grains, since the sphericity cannot be guaranteed for much finer fractions due to the production process.Of course one has to take care that the Knudsen regime prevails in the experiments, which requires a gas pressure below some Pa.
The generated packed samples were imported into SPARTA/DSMC as surface files of stl-format.The stl triangle meshes were constructed in MATLAB 3 from the output files (containing sphere centers) of the Baranau routines.Thus, a further preprocessing step was necessary to represent each sphere by surface segmentation of sufficient resolution (so as to avoid porosity errors).After import, the surfaces of spheres projecting beyond the simulation box are cut at the simulation boundary (cube surface).In contrast to SPARTA, the RW algorithm is based on ideal sphere geometry, so just the definition of the sphere centers and their diameter is needed as input.Either method (RW and DSMC) was used to obtain the height z the gas molecules reach when incident on the z = 0 face of the box.

Simulation results
As mentioned above, first RW as well as DSMC were used to determine the height particles reach in the sample.A large number of molecules (several millions) was used for this study step, incident to the bottom of the simulation box with Maxwellian velocity distribution.The fraction of molecules reaching height z is equal to the transmission probability T (z) of a layer of thickness z.The dependence of T on z obtained in this manner is exhibited in Fig. 6, where stars and circles mark RW and DSMC results, respectively.Though the two approaches are based on totally different algorithms, their results agree very well, confirming their accuracy.
In a second study step with DSMC, a constant incident flux of gas molecules is prescribed at z = 0, starting at time t = 0 with no gas in the sample.The applied gas species is N 2 , with a Maxwellian velocity distribution of incident molecules at a temperature of 296 K. Since the intermolecular collisions were switched off in the SPARTA script to force Knudsen diffusion, the applied pressure is irrelevant.Precisely speaking, the resulting density and flux distributions depend linearly on the pressure p 0 at the entrance (bottom side of the sample box), giving the same Knudsen diffusion coefficient independent of p 0 .The simulation was run until a steady-state flow was achieved within the sample.The evolution of the density and flux distribution in the sample with = 0.65 are plotted in Fig. 5.The density of the positive flux component, n + (z)/2, and the difference of forward and backward moving gas components, (n + (z) − n − (z))/2, are shown in addition to the total density n(z) = (n + (z)+n − (z))/2.Similarly, the flux  From the resulting gas density gradient dn/dz and flux J the Knudsen diffusion constant D K is calculated my means of (14).By applying (21) the half-transmission thickness L h is obtained, which finally is used in (9) to draw the transmission probability T (z).The corresponding curves are shown in Fig. 6 as solid lines (denoted as 'DSMC DK− >Lh' in the annotation).The results of this steady-state DSMC approach agree perfectly with the direct determination of the transmitted molecule fractions as obtained above with RW as well as DSMC (circle and star markers in Fig. 6).So the simulations provide a very good validation of the relation (21) between D K and L h .

Comparison of simulations with analytical models
In order to assess the accuracy of the analytical D K models, ( 9) is applied to the formulas ( 35)-(37), and the resulting functions T (z) are plotted in Fig. 6).So a comparison of the analytic models with the simulations is possible.
As already discussed by Asaeda et al. [4], their model overestimates D K significantly in a systematic manner (approximately by a constant factor of 1.4).Therefore they introduced the tortuosity correction factor q for compensation, just as usually done for the cylindrical Knudsen model (however, the two definitions τ 2 and q are different, since τ 2 contains the main part of the porosity dependence as discussed after (38)).
The relations (36) by Derjaguin and (35) derived in the previous section agree much better with the numerical results.On the average (35) performs best: At porosities above 0.6 the corresponding D K results deviate only about 3 %.The agreement deteriorates when approaching low porosity (high compaction), reaching about 10 % relative deviation at = 0.45.
The general Derjaguin model ( 36) has been calculated, as well, from the DSMC simulations.For this purpose, all chords (free paths λ between surface collisions) of the gas molecules are collected, for which the corresponding mean values λ and λ 2 are calculated and used in (36).This approach yields much better results than the Asaeda model (35).Since the latter can be derived from the former by approximating λ 2 = 2 λ2 , we can conclude that this equation is not well satisfied for the investigated packings.
Though the agreement of the models ( 36) and ( 35) with the numerical simulations is quite well, we have to state that in realistic packings the measured D K values are still lower, since in reality grains have angular surfaces, and even spherical glass beads have slightly uneven and rough surfaces.As a consequence, the specific surface of the pores is higher than the nominal value determined from the sphere diameters.So gas molecules undergo more collisions in the pores when they pass a layer of given porosity.This may but need not change the transmission probability through the layer.For instance tiny traps on the surface of very rough grains may cause a trapped gas molecule to collide many times before leaving the trap without changing the overall probability of passing the whole layer [32].The presented theoretical considerations in Sect. 2 are only applicable if the cross-sectional area porosity is sufficiently constant across the sample.However, this condition can be effectively violated when irregularly distributed tiny traps are present in a finite sample.In particular, when working with aggregates of finer grains, the pore space inside the aggregates can be regarded as traps, which may have only negligible influence on the transmission probability, although the gas molecules may stay (move) for a long time inside the aggregates.In order to quantify these effects better, the study of gas flow through hierarchical pore structures has attracted considerable attention recently [30,31,33].

Conclusion
The presented analytical approach facilitates the characterization of steady-state isothermal gas flow through plane layers of porous materials, when the distribution of pores is homogeneous on a scale greater than the typical pore size.It is assumed that the gas molecule reflection at pore surfaces is diffuse and no adsorption takes place.By considering the transmission probability T and reflection probability R = 1 − T of layer partitions it was shown that R/T is an additive quantity.In other words, its value equals the sum of the respective quantities R i /T i of a subdivision into n sublayers i = 1 . . .n along the net flux direction z.By taking the limit to infinitesimal sub-layers, it is recognized that L h = LT /R is independent of the layer thickness L and represents a property of the medium itself, which is revealed as the layer thickness for which half the incident gas molecules traverse the layer (the other half being reflected).
By applying a bi-hemispherical Maxwell distribution, splitting the gas density into the unidirectional components of forward and backward moving molecules, the link between gas density and flux components is made.This model was used by Asaeda et al. [4] to calculate the Knudsen diffusion constant D K of sphere packings.In the present article it is applied to relate T and L h to D K .This relation enables one to calculate the transmission properties of layers (given by T ) on the basis of D K models, or vice versa to find D K from measured fractions T of transmitted gas molecules.The relation proves useful in technical applications and laboratory experiments, but also in the evaluation of space science observations (comet outgassing), when one of the quantities can be determined.
By a further application of this stochastic approach to the reflections of gas molecules at the surface area elements in a thin sub-layer, an analytic model of D K is derived which differs from the former model by Asaeda et al. [4] by a constant factor 13/16.Numerical simulations with Random Walk and Direct Simulation Monte Carlo programs verify that the model derived here is more accurate than the former models.On the average, our layer model performs even better than Derjaguins general formula, which is based on the first two moments of the distribution of chord lengths (length between successive collisions of the gas molecules with the pore surface).The simulations further confirm the validity of the obtained relation between T , L h , and D K .In order that the presented probabilistic approach is valid, cross-sectional areas must be large with regard to grain or pore size, since the transmission probability must be independent of x and y (location on the plane normal to the net flux) on a macroscopic scale.More precisely, the probability definitions must be regarded as averages over all directions of incidence, and over areas much greater than the grain or pore size.
While the stochastic layer approach generally serves a better understanding of gas transport processes in porous media, it is of direct practical use, too.Since analytical formulas for the D K are often part of gas flow studies in the Knudsen regime, our models can contribute to a more accurate approximation of real scenarios in such cases [34,35].Similarly, the relation between L h and D K is a valuable tool, since it enables us to compare the results of simulations based on the former quantity with experiments and observations that are interpreted by the latter quantity.These merits are also motivation to generalize the models to more complex scenarios.For instance, an extension of the formalism to stratified media, where the medium properties (porosity and half-transmission thickness) change with height in the layer, is conceivable when the other mentioned preconditions are maintained.Corresponding work is planned in view of the many applications, e.g., laboratory measurements of granular samples and complex mesoporous materials, and gas transport in comet surface layers (which is not only relevant for the thermal evolution but also for dust entrainment of outgassing volatiles [36]).Of course, certain effects can scarcely be treated analytically and must be left to numerical simulations, e.g., irregular inhomogeneity or anisotropy of the pore  where the azimuth range 0 = max(0, ψ − π/2) < ψ i < min(ψ + π/2, π) = θ led to the limits of the inner integral in the second line.

Fig. 1
Fig. 1 Schematic of possible reflection and transmission scenarios of a gas molecule on its path through a two-layer porous medium with layer boundaries B 0 , B 1 , and B 2 .Different paths (a)-(e) are shown as examples

Fig. 2
Fig.2Molecule reflection at the surface of a sphere.The x axis is chosen in the plane spanned by the z-axis and the normal of the tangential plane (dotted) at the chosen collision point.The y-axis points away from the observer.The shaded area indicates the possible directions of incidence −e i of molecules with v z > 0 at the given point on the sphere.e i and the direction of reflection, e r , need not be in the shown cross section (xz-plane)

Fig. 5
Fig. 5 Density and flux evolution of the simulated gas flow through the sample of porosity 0.65 until attainment of steady-state conditions after about 20 ms

Fig. 6
Fig. 6 Transmission probability of a layer as a function of thickness z measured in grain diameters d g (equal to the fraction of molecules reaching height z from the entrance surface).The results of different models, as explained in the main text, are shown for three different porosities

Table 1
(21)-transmission thickness of homogeneous layers of packed spheres, as related to different models of D K via(21)

Table A1
Symbols for most important quantities and constants