An alternative procedure for modeling of Knudsen flow and surface diffusion

An alternative procedure for the calculation of impingement rate distribution and simultaneously the transmission probability in pores under Knudsen diffusion conditions is introduced. It is based on a combination of the finite difference method and a projection approach. Pore entrance and exit effects, and the influence of the pore length on diffusive fluxes are investigated. Later on, it is applied for a simultaneous Knudsen and surface flow system. In the model, the equation system is built without the independent flow and adsorption-desorption equilibrium assumptions. For the conditions investigated, the results indicate that if the surface flow rate is substantial, the independent flow and adsorption equilibrium assumptions become improper estimates for the behaviour of the system. The surface and gas flow rates, the impingement rate distribution and the surface coverage behave much more complex than the characteristics found with such assumptions.


Introduction
The common use of porous solids as adsorbents, membranes and catalyst-supports leads to the need of detailed understanding and modeling of transport in these systems.Achieving such an understanding not only involves the gas-phase diffusion but also the adsorptive properties of the material and the surface diffusion [1].Experimental studies have shown that surface diffusion contributes significantly to the total diffusive flux in both mesoporous and microporous systems [2,3].There are examples where the surface diffusion accounts for more than 50% of the total mass flow rate [4,5].
For the experimental determination of diffusion coefficients in such systems, both macroscopic and microscopic methods are available [1].For example, the most commonly used macroscopic technique is called the diffusion cell technique, which enables the user to measure surface fluxes of an adsorbed gas in a mixture while maintaining a very small surface concentration difference across the adsorbent in the cell [6].
Numerous studies have tried to elucidate the mechanism and the characteristics of surface diffusion [7].Today, there are basically three different approaches in modeling surface diffusion; mechanistic (hopping) model, random walk (Fick's law) model and hydrodynamic (slip) model [4,6].Fick's law approach is the most widespread in use [4,6] and assumes that the gaseous and the surface fluxes are independent of each other [6].Usually the connection between the gas phase concentration profile and the surface concentration profile is achieved via an adsorption equilibrium assumption [1,2,4,[8][9][10][11][12][13][14].
Although there are many studies, the diffusion in porous solids is still not clearly understood.For example: Hu et al. [15] report that the extracted surface diffusion is model-dependent and that the surface diffusion is complicated, Yang et al. [16] note that the Dusty-Gas-Model and various adsorption isotherm equations are not capable of representing transient behaviour of adsorbable gases and that the adsorbable gases exhibit a complex behaviour, Reyes et al. [1] remind that the differences frequently observed in diffusion parameters obtained by different experimental methods are not clearly understood and the discrepancies in diffusivity parameters have been discussed extensively in the literature.
The above complex behaviour and the discrepancies between different approaches have motivated the present authors to develop an alternative approach.It is claimed that the surface diffusion is significant only in the region where Knudsen diffusion prevails [17].Therefore as a starting point the authors investigate a system under purely Knudsen type gas diffusion.Secondly, the adsorption equilibrium assumption is commonly used without checking its validity.It was shown that physically it is not possible to have throughout equilibrium in such systems [18].Thus, it may be worthwhile to model the system without such an assumption to check its validity and to estimate errors involved in using it.
To be able to include all the above points in the modeling an alternative procedure is developed.The procedure is first used to reproduce pure Knudsen flow results to show its validity.Then, a system under simultaneous Knudsen and surface flow is taken into consideration.The model equations for this system is built without the independent flow and adsorption-desorption equilibrium assumptions.Finally, from the calculated results, the surface diffusion coefficient is back-calculated using these assumptions to rate the error in using them for such systems.

Theory for modeling Knudsen flow
When the mean free path of the gas is larger than the characteristic scale of the pore, the gas-surface interaction becomes a determining parameter in modeling of flow through pores.Knudsen argued that, under free molecular flow conditions, a cosine law of diffuse emission or reflection from the wall surfaces was the most reasonable assumption.He stated that each molecule is rejected with the same probability in any arbitrary azimuth, and the probability of a given angle of emergence is given by the cosine law.The direction in which a molecule is leaving the wall is completely independent of the incident one.[19], and [20] give good reviews of the Knudsen flow.It is also noted that the cosine law can be accepted as a correct basis for the evaluation of rarefied gas flow in both simple and complex vacuum systems as well as related fields [19].Nevertheless, it should be noted that under some conditions there are deviations from the expected cosine law [19], e.g., for structured surfaces deviations from Knudsen's law are shown to occur [21].
In the Cartesian coordinate system, the cosine law can be written as (see Fig. 1) where dn is the molecular flux through dω, N 0 is the total molecular flux from the surface element A, θ is the angle with the surface normal, dω(= sin θ dθ dφ) is the solid angle and φ is the azimuthal angle.Nevertheless, the implementation of the formula can sometimes be cumbersome and susceptible to mistakes [22].Additionally, the cosine law is not easily includable into the finite difference method in this form, but an alternative representation of this law presented below will allow for such an integration.Knudsen designed his experimental set-up such that there existed a surface element A on the inner surface of a spherical bulb from which the molecules have been scattered [23,24].These molecules are found to cover the inner surface of the spherical bulb homogeneously.The cosine law formula actually describes this distribution.In other words, if some molecules are scattered from a surface element A (which can be thought to be on a hypothetical sphere) (see Fig. 2) and directed to go through another piece of sphere surface, e.g., A 1 , then according to the cosine law, the ratio of these molecules to the total number of molecules scattered from A is equal to the ratio of the area of A 1 to the total area of the sphere.For example, for a cylindrical pore the fraction of molecules leaving an infinitesimally small surface element on the pore wall and passing through the cross-section at a distance h is of interest.The infinitesimal surface element A, the hypothetical sphere, and the cross-section at b are all shown for the case h = 2r in Fig. 3a.Projecting the pore cross-section (b, c) onto the sphere, center of projection being the surface element A, produces a cone, whose base is the cross-section (b, c), and whose tip is the surface element A (see Fig. 3b).Therefore, the sought area is the area on the sphere bounded by the intersection of this cone and the sphere.It should be noted that the ratio of the projection area to the total sphere area is independent of Per.Pol.Chem.Eng. the sphere radius, and thus for simplicity the pore radius and the sphere radius are taken be equal.To get an impression of the above mentioned situation, the sphere, the cone, and the projection area on the sphere can also be seen in Fig. 3b-e.
Since the molecules do not collide with each other and travel only by colliding with the pore wall, the calculation of fraction of molecules leaving a surface element in a particular direction is important.Considering that the fraction of such molecules can be calculated by using the projection of the cross-section, this approach can be named as projection approach.
There are basically two main quantities to be calculated for such a case.The fraction of molecules leaving a surface element and passing through a pore cross-section at distance h, F(h), and the fraction of molecules entering from the pore entrance and passing through the cross-section at distance h, G(h).The exact procedure using projection approach and mathematical formulae for the calculation of these functions are given in the appendix.If one discretizes the pore into many slices, and e.g.takes slice i into consideration, the system concerning the scattered molecules from a surface piece looks like as given in Fig. 4. In Fig. 4, the pore is divided into n slices with a constant thickness z, and the molecules leaving slice i in various directions are shown.Fi0 and Fi L are, respectively, the fraction of molecules leaving the pore from left and right ends.The f values seen in Fig. 4 represent the fraction of molecules leaving slice i and impinging on another slice.These f values can be calculated by taking the difference between two corresponding F values.In general, it can be defined as: and The f value between the i th and j th slice corresponds to The G(h) is defined as the fraction of molecules entering from the pore entrance and passing through the right boundary of the corresponding slice (see Fig. 5), that is, G(i • z) is the fraction of molecules entering from the left pore entrance and reaching the right boundary of the i th slice.Once again, the difference between the fraction of molecules reaching the right boundary and the left boundary of a slice gives the fraction of the molecules impinging on it, i.e. g(i) (see Fig. 5).Therefore, and also from Fig. 5, the fraction of molecules leaving the pore without impinging on it is A plot of F and G with respect to the normalized distance, h r por e , can be seen in Fig. 6.The total number of molecules entering the pore from the left entrance (Fig. 5) N lr z=0 , and from the right entrance N rl z=L , can be calculated from the kinetic theory of gases as  According to the kinetic theory, the number of molecules impinging on a unit surface is given by P• v 4RT , where R is the universal gas constant, T is the absolute temperature, and v = 8RT π M w is the mean average speed of gas molecules, with M w being the molecular weight of the gas molecules.Therefore when multiplied with the pore cross-section, πr 2 por e , it gives the amount of molecules entering from the pore entrance.

Modeling
A general model including both surface and Knudsen flow will be set-up first.Since sole Knudsen flow will be a simplified version of this general model, its equations can be derived by simplifying these more general equations.
Mostly, when a porous system is modeled, the facial outer surface area is not included in the model.Although for a purely gas phase diffusion system not being important, for the case where there is surface diffusion, such a surface area becomes important to determine the boundary conditions for the surface flow.Such a model is given in the work of Argonul et al. [18] and is going to be used in this work too.
A pictorial representation of the system, a cylindrical pore with both gas phase and surface flow, is given in Fig. 7.In the figure , C A is the gas phase concentration of component A, G A is the surface concentration, r ads and r des stand for adsorption and desorption rates and the F Sur f and F gas respresent the surface and the gas phase flow rates, respectively.In Fig. 7, the facial outer surface area surrounding the pore entrance (see Eq. 33 below) is the area on the left-end and right-end of the solid substance facing the bulk of the gas.Since F Sur f 0 originates from the left outer surface (los) and F Sur f L terminates at the right outer surface (ros), the boundary conditions for the surface flow are determined through the mass balances at the outer surfaces.
The system is modeled under the following conditions: • Steady-state flow • No collisions between the molecules, i.e., pure Knudsen diffusion in gas phase • The flux ([ mol/(ar ea • time)]) of the incoming molecules at the pore entrance is homogeneously distributed over the entrance cross-sectional area, the velocities correspond to the average of the Maxwell distribution corresponding to a given temperature and the flow direction of the molecules follow the cosine law.
• Diffuse scattering (with cosine law) for the collisions with the walls • Constant surface diffusion coefficient • Langmuir type adsorption • Monolayer surface diffusion, i.e., surface diffusion under the condition that the surface is covered below the monolayer adsorption capacity [4].
Under these conditions, total number of molecules impinging on slice i, N imp (i), can come from three different sources: molecules from the left pore entrance, N lr z=0 • gi ; molecules from right pore entrance, N rl z=L • gn−i+1 ; and molecules from other slices, ) It should be noted that since a constant slice thickness is used, i.e. z = const., the g values for the left entrance can be used also for the right pore entrance by making simple index switching.Since the i th slice from the left is the (n − i + 1) th slice from right, gn−i+1 should be used to calculate the fraction of molecules impinging on i th slice that are entering from right pore entrance.The Eq. 7 is only for the i th slice, if all the equations for n slices are written, they can be combined to give a single equation in matrix form (note that g f li pped (i) = g(n − i + 1)): where . . .
) and fn×n is a symmetric Toeplitz matrix.In general, a Toeplitz matrix is only perm-symmetric [25].
A mass balance for slice i in the pore can be visualized as in Fig. 8, where R ads,i and R des,i are the adsorption and desorption rates [mol/s], and N dir,i represents the amount of molecules reflected directly [mol/s], i.e., without adsorbing.
From Fig. 8, it can be seen that the total amount of molecules leaving the slice, i.e., scattered, is equal to the sum of molecules desorbed and molecules directly reflected: additionally, the impinging molecules are either directly reflected or are adsorbed.
Direct reflection of molecules can be due to three reasons; first, molecules impinging on a surface point that is not an adsorption site, second, molecules impinging on a surface piece that is appropriate for adsorption but already occupied by another adsorbed molecule and third, if the sticking coefficient, s coe f , is not unity, molecules impinging on free adsorption sites that are not adsorbed.If β ads denotes the ratio of the surface area capable of adsorption (area of the adsorption sites) to the total surface area (obtained, for example, from BET measurements), and the adsorption is of Langmuir type, then the rate of directly reflected molecules, N dir , can be related to the above three points, respectively, as (θ i ≡ surface coverage): In a simplified form Eq. 13 becomes: The β ads can be estimated from the following formula where σ A is the area of one adsorption site, which can be estimated to be equal to the area of the adsorbing species [m 2 /molecule], q mono is the monolayer coverage [mol/g cat ], N Av is Avogadro's number [molecules/mole], S B E T is the B.E.T. surface area of the porous structure [m 2 /g cat ] and G total is the surface concentration of appropriate adsorption sites [mol/m 2 ].
According to Langmuir type isotherm, the desorption rate is proportional to the desorption rate constant, k des [mol/s/m 2 ], and to the surface coverage.
where k des , [mol/s], is the desorption rate constant based on the surface area of one slice.Thus, if Eq. 14 and 16 are inserted into Eq.11 and written in matrix form for the whole system, one ends up with: where " ×" represents element by element multiplication (Hadamard product).
Recalling the definition of N scat , i.e., Eq. 11, and making a mass balance for the surface element in Fig. 8, one can write that: where the F Sur f represent the surface molar flow rates.Thus Eq. 18 can be written for the i th slice, For surface diffusion, Fick's first law combined with the finite difference approach reads as follows: From Fig. 9 it can more clearly be seen that the F Sur f 0 and F Sur f n are dependent on the left and right outer surface (los, Fig. 9.A plain representation of the surface flow in the porous substance that is divided into n slices.ros) coverages.If one takes the distance between los and the first slice, and between the last slice and ros as half of the slice thickness, i.e., z/2, then the corresponding surface flow rates become: The mass balance for the first and the last slice then turn out to be: In general, it can be written that Consequently, Eq. 8 and 25 can be combined to give which gives the dependence of N scat on surface coverages (hidden in F Sur f ).
The gas phase flow rate inside the pore can be calculated based on the projection approach as follows (note that In words that means: the molecules reaching the crosssection (at position i • z) from the left entrance (N lr z=0 • G(i • z)) minus the molecules coming from right entrance (N rl z=L • G ((n − i) • z)) plus the molecules coming from slices on the left of the cross-section ( F(1 minus the molecules coming from slices on the right ( F(1 : gives the net gas flow through that cross-section.
At steady state, the conservation of mass principle requires that the total flow rate should be constant throughout the whole pore.This requires the validity of the following equation at any cross-section: Using the total flow rate at the left entrance (F total 0 ) as basis, one can set up n independent equations describing the principle of constant total flow (Eq.28) Eq. 17 can be rewritten as follows: The above two equations (29,30) form the basis for the solution of the model.In addition to these, another equation can be set up which accounts for the fact that at steady-state there can be no accumulation of mass on the surface.This means that the net adsorption (i.e., the difference between adsorption rate and the desoption rate) over the whole surface should be zero.That covers both the outer facial surfaces and surface inside the pore.It should be noted here that an overall adsorption/desorption equilibrium satisfies this condition but is a special case and is shown to be generally not possible [18].Thus, If the surface concentration on the los is assumed to be homogeneous, i.e., no concentration gradient on los, owing to mass balance, the rate of net adsorption on that surface would be equal to the surface flow rate (see Fig. 7): The A f ace solid/ por e is the facial outer solid surface area per pore and can be calculated by making use of the porosity , of the porous structure and the pore cross-sectional area as: The rate of adsorption and desorption at los according to Langmuir approach are: The adsorption rate is equal to the number of molecules impinging on the surface multiplied by the fraction of surface available for adsorption, times the fraction of free adsorption sites, times the sticking coefficient.Equating this to the known Langmuir adsorption rate would lead to the adsorption coefficient.
and consequently If one has the adsorption equilibrium constant K ads , then the calculation of k des is straight forward

The model equations for sole Knudsen flow
When there exists no surface flow, at steady-state, scattering rate must be equal to impingement rate.From Eq. 25 it leads that additionally one can assume that there exists absolute vacuum on the right-hand side of the pore, then the main equation to be solved (Eq.26) becomes An important quantity for Knudsen flow is the transmission probability.The transmission probability w is equal to the fraction of molecules leaving from the right pore-entrance, i.e.N out right , divided by the total number of molecules entering from the left.
The N out right can be calculated by summing up the molecules passing through the pore without impinging on it (see Fig. 5) and the molecules scattered from the pore slices in the direction of the right pore-entrance (see Fig. 4): where Fi L 1×n is a row vector with elements corresponding to the fraction of molecules scattered from slices in the direction of the right pore-entrance.

Summary of the general model
A single pore (or a simple parallel pore structure) with Knudsen flow accompanied by simultaneous surface flow with Langmuir type adsorption is modelled at steady-state.It should be noted that the generated model does not make the assumption of adsorption-desorption equilibrium, and it incorporates automatically the variation of flow rates with pore length and the so called entrance and exit effects.Any impinging molecule is either adsorbed or directly reflected.Scattered molecules are taken to be the sum of desorbed molecules and directly reflected molecules.The scattering is assumed to follow the cosine law of reflection.For surface flow, Fick's law is assumed to be valid with a constant surface diffusion coefficient.Finite difference approach is followed and the pore length is discretized into n slices.Through mass balances, the necessary equations are set up.
The system of equations to be solved is Eq.29 and 30 with surface coverage (θ) as unknown.Required auxiliary equations are Eqs.8, 20-22, 26, 27, 33-37 and 39-40.The solution should also satisfy Eq.31, therefore this equation can be used to check the validity of the obtained solution.Details of the solution procedure are given below.

Physical constants used for the general model
A sample set of constants is taken from Chen and Yang [26] where the values have been extracted from the experimental measurements of Gilliand et al. [27].As an example, the system propylene (C3H6) in (Vycor-) glass is taken, which was found to exhibit monolayer adsorption, which corresponds to our model.The data is given in Table 1, and other parameters used are shown in Table 2. [−] [m] * Poling et al. [28] 7 Other models used for comparison Two other models are used for comparison with the model in this work, and the models are numbered with increasing detail.The first one, i.e., model I, assumes independent gas phase and surface flows with established adsorption equilibrium, and additionally linear impingement rate distribution inside the pore between the reservoirs with different pressures at the two ends.
The second model, i.e., model II, also assumes independent flows and adsorption equilibrium, but uses the projection approach for pure Knudsen flow to calculate the impingement rate distribution inside the pore.The model set up in this work is labelled as model III, and as explained in the previous sections it does not make any of the above assumptions.The gas phase and surface flows are related to each other through adsorption and desorption rates that are separately calculated.

Solution of equations
The equation system of model III is normalized by the theoretical gas phase flow rate calculated by the traditional Knudsen diffusion coefficient (i.e., nor mali zer = − 2 3 r por e L por e • (P right − P le f t ) • v RT • πr 2 por e ).The normalized equation system is then solved using MATLAB (7.3.0,R2006b).The function 'fsolve' with 'LargeScale' function set to 'on' is used to solve the system of nonlinear equations.The equilibrium distribution of model II is used as an initial guess for the solution.To check the convergence behaviour of the model, a sample system was simulated by modified initial guesses (which were taken to be very different from each other) and the system converged to the same solution for various runs made.The calculations for the sole Knudsen flow case is also done by MATLAB.The flow diagrams for the two different cases, i.e. sole Knudsen flow and Knudsen flow with surface flow are given in Figs. 10 and 11.
A note should be made here concerning the used slice thickness.The pore is originally assumed to be discretized into smaller elements than the slices, which are named sub-slices.And the F and G values for them are calculated before hand and saved.These finer division of the pore provides higher resolution for the interaction of the pore slices within each other and with the pore entrances.But when it comes to the calculation of, e.g.N imp and surface coverages, the equation system becomes unnecessarily large if one uses the sub-slices.One does not need such a high resolution for the N imp and surface coverage values, which are practically constant for a series of sub-slices.Consequently, the sub-slices can be bundled together into what is now called slices.For example if the relative slice thickness is chosen to be 0.01 ( z = z r por e = 0.01), and the relative subslice thickness is 0.002 (δ z = δz r por e = 0.002), five sub-slices are to be bundled together.The F and G values for the slice can then be calculated by taking the average of these values for the sub-slices.The F and G values calculated in that way are more precise than the calculation using the corresponding slice, and also, for example, the matrix fn×n is 25 times smaller than the case of using sub-slices.One just needs to choose appropriate values for the z and δ z considering the system at hand.

Sole Knudsen flow case
As an example a pore with a length to radius ratio of ten is simulated with various slice thicknesses (with a relative sub- slice thickness, δ z = δz r por e = 2 • 10 −3 ).The results are tabulated in Table 3.It can be seen that as the relative slice thickness z = z r por e gets smaller, the transmission probability becomes more precise.This also means that one can choose z (and δ z) accordingly to reach the desired level of accuracy.
To check the accuracy of our program and the model, calculated values from the projection approach are compared with values from the literature obtained by various methods (see Table 4) .Column one shows the values calculated by the Knudsen formula 1 , column two presents the values of DeMarcus taken from [29], column three are values of Berman taken from [30] and column four are simulation results reported by [31].The last two columns are the values calculated in this work.Two different slice thicknesses are used for the calcula-1 Although Knudsen formula is valid for very long pores, it is included to show the extent of error made when it is to be used as an approximation.
Tab. 3. The effect of slice thickness on the transmission probability for a pore of length of ten times the pore radius with sub-slice thickness of 2 tions in order to be able to estimate the last significant digit.A very good agreement of our results with the values given in the literature can be ascertained.
For the impingement rate results, the table from [32] is used for comparison (see Table 5).They have compared the impingement rate at the entrance pore surface found by analytical approximations in the literature (first four columns) with their numerical solution (fifth column).The values using projection approach are found by fitting the calculated impingement rate values to a cubic spline curve and then extrapolating this curve to the pore entrance.
All the above results indicate that the projection approach delivers accurate results.Here, it should be noted that the method and its principles are not only applicable to cylindrical pores but to any cross-sectional shape.The cylindrical pores is chosen as an example due to its symmetry, which highly simplifies the calculations, and its common application in modeling diffusion.

Surface flow and Knudsen flow, general model
A pore length of L por e r por e = 20 has been chosen as an example.The pore is divided into slices of relative thickness of z = 0.05.The pressure on the left-hand side of the pore was taken to be higher than the pressure on the right-hand side by 20 kPa, i.e., a constant pressure difference across the pore length, for all the runs.
Total flow rate (Knudsen plus surface flow rate) and the flow enhancement with respect to the case of purely Knudsen flow are tabulated in Table 6.
Tab. 6.Total flow rate and flow enhancement with respect to pure Knudsen flow case for various pressures It should be noted that the flow enhancement is neither linearly dependent on the average pressure nor constant like the Knudsen flow rate.
The two extreme pressures can be analysed as examples, since the behaviour of the system in between can be deduced from these two extremes.
For the case P le f t = 30k Pa the impingement/scattering rate distribution is shown in Fig. 12a.The broken line with dots represents the linear pressure distribution inside the pore, i.e., model I.It starts from the left-hand side bulk pressure value and decreases linearly to the right-hand side bulk pressure value.The solid curve represents the pure flow case and consequently the independent flow case (i.e., model II).The broken curve and dotted curve are the results of model III and represent the impingement rate and the scattering rate, respectively.Apparently, the behaviour of this particular system is different from the assumption of independent gas and surface flows in equilibrium.Although the calculated scattering and impingement rates seem to be close to each other, i.e., around equilibrium, model III allows the curves to shift away from the expected equilibrium curve, i.e., model II.It can be noted that although the model I and model II give symmetric curves, for model III the resulting distributions not symmetric around the mid-point of the pore anymore.This shift and the unbalance are the results of the effect of the surface flow on the gas phase flow.
The relative difference between the adsorption and desorption rates can be better seen in Fig. 12b, where the adsorption/desorption rate difference (ADRD) is plotted, normalized by the adsorption rate versus pore length.
The triangles in the figure represent the values on the left and right outer surfaces, where actually the biggest difference between adsorption and desorption is observed.It can be noticed that at los around 8.7% and at ros around 24.5% difference is created.This indicates that there is a considerable transfer of species between the gas-phase and the surface-phase at these regions.Inside the pore, the entrance and exit regions have nonequilibrium conditions, but the mid-region can be said to be at quasi-equilibrium.But as can be seen in Fig. 12a, this quasiequilibrium values cannot be calculated by simply assuming independent flows with equilibrium distribution.The system behaves close to equilibrium but far from the state found with independent gas flow and equilibrium assumption combination, and also behaves differently.
Nevertheless, the behaviour in Fig. 12b was expected [18], the system is anticipated to have a net adsorption region first, then a local equilibrium point after that a net desorption region.Besides, the los was expected to have net adsorption and the ros net desorption.
The adsorption and desorption rates are continuous inside the pore, but there needs to be a jump between the values just at the entrance of the pore and the outer surfaces due to the jump of the impingement rate (see Fig. 12(a), the model III curves do not start from 1 and end at 1/3) at the same region.This is the reason for the discontinuity at Fig. 12b at two ends of the pore, between the curve inside the pore and the two end points at the outer surfaces.Note also the asymmetry at the pore ends (Fig. 12b) due to the different pressures in the reservoirs.
The surface coverage profiles for all three models are given in Fig. 13.Unlike the adsorption and desorption rates, the surface coverages and also the surface flow rate (Fig. 14a) are not discontinuous at the two ends.As can plainly be seen in the Fig. 13, neither model I, nor model II are good estimates of the behaviour of this particular system.Fig. 14 is maybe the most interesting graph for the system.
As presented in the Fig. 14a, the surface flow rate originates from the left outer facial surface and then increases along the pore length.This increase is most pronounced at the entrance region of the pore.In the deeper parts of the pore, the increase is small but still existent.These effects are due to net adsorption along this pore portion.Since the total flow rate is constant (steady-state flow), the gas phase flow rate decreases accordingly (see Fig. 14b).Close to the end of the pore, the surface starts to desorb more molecules than it adsorbs, and the surface flow rate begins to decrease after some points in the pore, and finally comes to the right outer surface value.The amount of surface flow reaching the ros is given back to the bulk of the gas by net desorption from the outer surface, which then causes the big difference between adsorption and desorption there.
It should be noted here that, for example for the other extreme case, i.e., P le f t = 150k Pa, the impingement rate and surface coverage are very close to the equilibrium distribution (see Fig. 15a-b).But the surface flow follows a similar behaviour (Fig. 15c) as in the previous case.Here model II and model III behave similarly except for the entrance and exit regions of the pore.A notice should be made here concerning the surface flow rates.In the latter case the surface flow enhancement is only 13.17%, but for the first case (P le f t = 30k P A) it was 95.74%.That is due to the shape of the adsorption isotherm.Since the Langmuir isotherm is not linear, a constant pressure difference does not lead to a constant surface concentration gradient, thus   Since the surface flow is low, the gas phase flow does not change much in magnitude along the pore (see Fig. 15d) even if it follows a similar behaviour as in the previous case.Consequently, it may be said that if there is considerable surface flow, at least for the pore system investigated, the independent flow with adsorption equilibrium assumption is not representative anymore.

Back calculation of the surface diffusion coefficient
It may be interesting to calculate the surface diffusion coefficient from the simulation results, assuming independent flows.Since independent flow with equilibrium assumption is commonly used in the literature, it may be useful to have an idea of the error involved employing it for such a system.The following common assumptions are made for such a calculation: 1 pure Knudsen flow in gas phase 2 surface flow and gas phase flow are independent and thus can be calculated independently 3 the pressure difference across the pore is small (20kPa) 4 the gas phase and the surface are in equilibrium The flows are then simply where gas phase flow is constant due to constant pressure difference and has a value of 2.4678 • 10 −15 [mol/s] (based on model II) and the surface flow rate is defined by Fick's first law of diffusion as: Using the Langmuir isotherm, the dependence of surface concentration gradient on pressure gradient can then be calculated as: since the pressure difference across the pore is small, the gradient can be replaced by the difference, and an average pressure can be used P avg = P le f t +P right 2 in the denominator.Consequently, the surface diffusion coefficient can then be estimated by combining and rearranging Eqs.46 and 47: The relative value of the estimated surface diffusion coefficient in this way against the average pressure is plotted in Fig. 16.
There are two important points that Fig. 16 leads to.First, the assumption of independent flows and adsorption equilibrium in this particular case leads to errors roughly between 15% and 25% in the surface diffusion coefficient.Second, the estimated surface diffusion coefficient values follow a definite path with respect to average pressure and may mislead to the rash conclusion that the surface diffusion coefficient is not constant but concentration dependent.It is already known that various experimental papers about surface diffusion conclude similarly, i.e., surface diffusion coefficients are concentration dependent.Now the question arises if such a conclusion is correct or it is merely a misconception being the result of the assumptions made.It should be noted that the system may behave differently for different values of parameters, such as surface diffusion coefficient, adsorption constant, monolayer loading, etc.. Thus, the results and discussion should not be taken to be describing the general behaviour of all such systems.Therefore, it is not possible to give just one answer for such a question, instead one has to refer to the particular conditions.

Conclusions
An alternative procedure for the calculation of impingement rate distributions and simultaneously the transmission probabilities under Knudsen conditions in a cylindrical pore has been introduced and shown to give accurate results.In this approach the fraction of molecules leaving a surface element and passing through a cross-section of the pore at a distance h is found by taking the projection of the cross-section onto a hypothetical sphere generated from a defined surface element.This is basically another way of interpreting the well-known cosine law of scattering.
The procedure is then used for modeling a system under simultaneous Knudsen and surface flow.The model also makes use of the outer surface area of the solid to calculate the boundary conditions for the surface flow and is free of the adsorptiondesorption equilibrium assumption.
The results for the investigated conditions indicate that if the surface flow rate is significant then the independent flow and adsorption equilibrium assumptions are not able to describe the behaviour of the system adequately.The surface and gas flow rates, the impingement rate and the surface coverage distributions are considerably different from the prediction made based on such assumptions.For the case, where the surface flow rate is highest, a 25% difference was found between the adsorption and desorption rates on the right outer surface, indicating the non-equilibrium condition and the importance of the usage of the outer surface in modeling.
Additionally, the calculated surface diffusion coefficients using the independent flow and adsorption equilibrium assumptions from the results were underestimated by 15% to 25%.Moreover, these estimated diffusion coefficients are seemingly concentration dependent, although a constant surface diffusion coefficient was used.This arouses questions about (macroscopic) experimental studies resulting in concentration dependent surface diffusion coefficients.In such works, similar assumptions, i.e., independent flow and adsorption equilibrium, are commonly used, but the importance of the neglected effects is a priori unknown.
Anyhow, one cannot conclude from the present calculations that for every system these assumptions will lead to erroneous results.Many systems and various conditions should be investigated to give more concrete estimations of error and the behaviour of the representative systems.

F
A sketch of the coordinate system, sphere and cone are given in Fig. 17a.
The equations describing, respectively, the sphere and the cone are: The intersection curve between the sphere and the cone can then be found by equating the above two equations.The resulting formula for z is then: The above equation is defined only for 0 x int x max , and the y int equation (parameterized with respect to x) can be found by inserting the above equation into Eq.49.The value for x max can be found by using similar triangles in Fig. 17a as: after inserting Eq.51 and rearranging, one gets: If one takes a slice at an x int , as seen in Fig. 17a, one gets an arc on the sphere on the y-z plane (Fig. 17c).The r i in the figure is the radius of the circle cut from the sphere at x int .It is defined as (from the equation for the sphere, i.e.Eq.49): The cosine of the angle α that subtends the half of the arc is then: If one recalls Archimedes' rule that the surface area of a slice of a sphere and of a cylinder (whose radii are the same) are the same; instead of calculating the surface area on the sphere, one can calculate the surface area on the cylinder.The surface area on the cylinder can be calculated by adding up (i.e., integrating) the lengths of the arcs on the cylinder.The length of an arc is the radius times the angle it subtends.Since the angle calculated above (for the sphere) is also the same for the cylinder, it follows then that the length of the arc on the cylinder is: Then the integral required to find the surface area is straight forwardly: A special case would appear if the distance between the base of the cone and the sphere center (i.e., h) is less than twice the radius of the cone (i.e.r c ).In that case, a cap of the sphere comes additionally to the above calculated area as being part of the sought area (see Fig. 17b).The area of a cap of a sphere can easily be calculated if the height is known.In this case, the height of the cap is: and then the area of the cap: Thus the total area enclosed on the surface of the sphere can be generalized as: finally B Calculation of projection area for G The molecules enter the pore from the left pore entrance and reach the cross-section at a distance h (see Fig. 18a).The vertical line at this distance is the base of the cylinder.
The equations describing the sphere and the cone are, respectively: For the determination of the cone equation, the center line equation and the radius of the cross-section of the cone at an arbitrary x (found by using similar triangles) is used.The intersection curve between the sphere and the cone (found by equating Eq.62 and 63) is: It should be noted that Eq.66 is only valid for x int 0 and x min ≤ x int ≤ x max .In order to find x min and x max , one needs the equations for the top-line and for the bottom-line, which are: Inserting these two equations into the sphere equation, one finds the intersection points of the sphere and the two lines.The xcoordinates of the intersection points are zero and: x min = 2r s From Fig. 18b, it can be seen that the cosine of the angle α is: Note: The minus sign in the above equation is due to its position with respect to the origin used to define the angle.The origin for the angle is (x int , 0, r h ) and the negative z direction is defined to be the positive direction for the corresponding axis.Since in this case the angle turned out to be greater than π/2, the cosine of this angle is supposed to be a negative value.Therefore, the length (z int − r h ) has to be multiplied with minus one.The length of the arc that would be built on a surface of a similar cylinder would then be (see Appendix A -Calculation of Projection Area for F): L ar c = 2r s • ar ccos r h − z int r i (72) The radius of the arc can be found by utilizing the equation of the sphere (Eq.62): Consequently, the area of the piece of the sphere bounded by the intersection curve between x min and x max is then: and if it exists, the area of the cap left behind to the right of the x max is: Thus the total area enclosed on the surface of the sphere for this case also is: The fraction of this area to the total sphere area can then be easily calculated: Gh is the fraction only at the circle with radius r h at the entrance.The overall fraction can then be calculated by the following formula: Modeling of knudsen flow and surface diffusion 53

Fig. 1 .
Fig. 1.Representation of cosine law.A is a plane surface element, θ represents the angle with the surface normal, dω is the solid angle.The molecular flux through dω, i.e., dn, can be related to N 0 , the total molecular flux from the surface piece A, as follows: dn = N 0 π cos θ dω.

Fig. 2 .
Fig. 2. Alternative representation of the cosine law.The fraction of molecules scattered from A and passing through A 1 is found by the ratio of the area of A 1 to the total sphere area, dn N 0 = A 1 4πr 2 spher e

Fig. 3 .
Fig. 3.The surface element A, the imaginary sphere (radius r ), the crosssection at h (h = 2 • r ) are shown in a pore of radius r .(a) A 2-D representation of the system (b) A 3-D representation of the system (c) A 3-D representation of the projection area on sphere (d) Side-view of the projection area (e) Top-view of the projection area.

Fig. 6 .
Fig. 6.The function F and G against the normalized distance between the emission point and the cross-section of interest.The distance is normalized with respect to the pore radius, and F(0) = 0.5 and G(0) = 1.

Fig. 4 .
Fig. 4. The distribution of flow from a slice to the other slices and to the two ends of the pore.

Fig. 5 .
Fig. 5.The distribution of the flow from the (left) pore entrance between slices and pore exit.

Fig. 7 .
Fig. 7.An overview of the pore and the outer surfaces with the transport processes involved.Outer surfaces are the solid surfaces facing the bulk of the gas on both sides of the figure (G A,0 and G A,L are the surface concentrations on these surfaces).

Fig. 8 .
Fig. 8.The pictorial representation of the mass balance for a control volume at the surface of the pore.

Fig. 11 .Fig. 10 .
Fig. 11.The flow diagram of the computer program for Knudsen flow with surface flow (model III).

Fig. 12 .
Fig. 12.(a) The impingement and scattering rate distribution inside the pore.All the are normalized by the impingement rate of the gas at a pressure of P le f t .(b) The percentage difference between adsorption and desorption rates

Fig. 15 .
Fig. 15.(a) The impingement and scattering rate distribution (normalized by the impingement rate of the gas at a pressure of P le f t ) inside the pore.(b) Surface coverage profile.(c) Surface flow rate profile.(d) Gas phase flow rate profile.(for the case P le f t = 150k Pa)

Fig. 16 .
Fig. 16.Plot of estimated surface diffusion coefficient normalized by the actually used surface diffusion coefficient against the average pressure inside the pore.The estimation involves the independent flow and adsorption equilibrium assumptions.

Fig. 17 .
Fig. 17.Auxiliary figures for the calculation of the projection area in F calculation.(a) Sphere, cone and the origin of the coordinate system used, (b) the cap built on the right of the x max if h < 2r c , (c) the arc built on the sphere at the plane x int .

Fig. 18 .
Fig. 18.Auxiliary figures for the calculation of the projection area in G calculation.(a) The sphere tangent to the entrance at position r h and its relation to the cone whose base is at a distance h from the entrance, (b) a sample arc built between x min and x max .
NomenclatureA: in subscript denotes a property of species A A : denotes a small surface element A f ace solid/ por e : outer surface area available per pore, [m 2 /pore] AD R D : abrr.for Adsorption Desorption Rate Difference C A : concentration of A in gas phase, [mol/m 3 ] D Sur f A : surface diffusion coefficient, [m 2 /s] D Sur f est.: back-calculated surface diffusion coefficient using the adsorption-desorption equilibrium assumption, [m 2 /s] f : function giving the fraction of molecules scattered from pore surface and impinging onto a pore slice at a distance F : function giving the fraction of molecules scattered from pore surface and reaching pore crosssection at a distance F f li pped : the flipped version of the vector F, i.e., F f li pped (i) = F(n − i + 1) F gas : gas phase flow rate, [mol/s] F Sur f : surface flow rate, [mol/s] F total : total flow rate, i.e., F gas + F Sur f , [mol/s] g : function giving the fraction of molecules entering from the pore entrance and impinging onto a pore slice at a distance G : function giving the fraction of molecules entering from the pore entrance and reaching pore cross-section at a distance g f li pped : the flipped version of the vector g, i.e., g f li pped (i) = g(n − i + 1) G total : total concentration of adsorption sites on the surface, [mol/m 2 ] G A : concentration of A on the surface, [mol/m 2 ] I n×n : identity matrix of size n k ads : adsorption rate constant, [mol/s/m 2 /Pa] K ads : adsorption equilibrium constant, [1/Pa] k des : desorption rate constant, [mol/s/m 2 ] k des : desorption rate constant for a slice, [mol/s/slice] los : abbr.for left outer surface, in subscript it denotes the value on the los L por e : length of pore, [m] M w : molecular weight, [kg/mol] n : total number of slices N lr z=0 : rate of molecules entering from the left pore entrance, [mol/s] N rl z=L : rate of molecules entering from the right pore entrance, [mol/s] Per.Pol.Chem.Eng.
Comparison between calculated transmission probabilities w (for δ z = 2 • 10 −3 ) with various approaches from literature.