- Email: [email protected]
Simulation of gas diffusion in highly porous nanostructures by direct simulation Monte Carlo
Chemical Engineering Science 105 (2014) 69–76
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces
Simulation of gas diffusion in highly porous nanostructures by direct simulation Monte Carlo Jochen A.H. Dreyer a,b, Norbert Riefler b, Georg R. Pesch b, Mirza Karamehmedović b, Udo Fritsching b, Wey Yang Teoh a, Lutz Mädler b,n a b
Clean Energy and Nanotechnology (CLEAN) Laboratory, School of Energy and Environment, City University of Hong Kong, Hong Kong SAR Foundation Institute of Materials Science (IWT), Department of Production Engineering, University of Bremen, Germany
H I G H L I G H T S
G R A P H I C A L
A B S T R A C T
Direct Simulation Monte Carlo of gas diffusion with the open source solver dsmcFOAM. Improved accuracy by implementation of the variable soft sphere model. Correction of the particle flux accumulator to avoid pulsed gas initialization. CO diffusion in porous anisotropic nanostructures as formed with FSP. Base for describing adsorption, desorption and chemical reactions.
art ic l e i nf o
a b s t r a c t
Article history: Received 4 September 2013 Received in revised form 17 October 2013 Accepted 22 October 2013 Available online 6 November 2013
A Direct Simulation Monte Carlo (DSMC) method is utilized to simulate gas diffusion in nanoscaled highly porous layers. An open source solver has been extended with the variable soft sphere (VSS) binary collision model and the inflow boundary model was adjusted for small numbers of DSMC particle initialization. Comparison with the analytical diffusion equation illustrate the improvement of the VSS model compared to the variable hard sphere model (VHS). Subsequently, several highly porous particle layers (gas sensors synthesized by flame spray pyrolysis and isotropic layers) build up by 10 nm particles have been investigated. Results for DSMC gas diffusion in the porous structures are in agreement with the well established dusty gas model (DGM). However, while DGM requires measurements or estimations of pore sizes, porosity, and tortuosity and furthermore is limited to homogenous layers, the present contribution shows significant advantages of DSMC in describing gas diffusion in nonisotropic porous structures. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Nanoparticle aggregate films Gas diffusion DSMC Dusty gas model High porosity Gas sensor
1. Introduction Diffusion of gas molecules into porous structures is a phenomenon which occurs in many environmental or technical systems such as gas sensors and catalyst beds. The understanding of macroscopic effects in such systems often requires modeling at the mesoscopic scale. To describe rarefied gases or systems with
n
Corresponding author. Tel.: þ 49 421 218 51200; fax: þ 49 421 218 51211. E-mail address: [email protected] (L. Mädler).
0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.10.038
small pores (Knudsen number Kn 4 0:1), e.g. mesoporous systems (pore-diameter between 2…50 nm), kinetic theory has to be applied through solving the Boltzmann equation. Early approximations of this equation for pure diffusion without walls led (besides the Boltzmann H-theorem) to the Chapman–Enskog perturbation method, followed by the Grad–Zhdanov moment technique. The linearization of the Boltzmann equation was introduced later by Bhatnagar, Gross and Krook (Cunningham and Williams (1980)) and serves as the basis of many lattice Boltzmann model (LBM) simulation methods (Succi, 2001). In contrast, the dusty gas model approach (DGM) uses the
70
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
non-linear Chapman–Enskog or the Grad–Zhdanov constitutive equations of diffusion, where wall effects are considered by collisions of the gas molecules with stationary rigid spheres (Cunningham and Williams, 1980). Another approach to calculate diffusive and other kinetic properties is the direct tracking, i.e. calculation of tracks and collisions of gas molecules. Molecular dynamics (MD) is an often applied method where each atom or molecule is individually tracked, however, with the expense of very high computational costs if larger problems are considered. A good compromise between the required computing power and problem size is the Direct Simulation Monte Carlo method (DSMC) (Bird, 1994). In this method, the molecules of a volume cell are united into different parcels. The translation and the collision process are decoupled in every time step. In the collision process, only a part of the molecules are chosen according to their collision cross section. These collision partners represent statistically the ensemble of all molecules within the simulation volume. This work describes the implementation of the variable-softsphere (VSS) collision model (Koura and Matsumoto, 1991) in an existing DSMC solver of OpenFOAM. The results are compared with analytical solutions of the diffusion equation. Beyond that, the diffusion processes into gas sensor layers are investigated based on this implementation and compared to the well established DGM method.
2. Theoretical models
size range of the pore diameter. This case occurs for Knudsen numbers Kn Z 0:01 with Kn ¼ λ=daver: , where daver: is the average pore diameter and λ the mean free path. In case of non-adsorbing gas molecules (Petropoulos and Papadokostaki, 2012) and for Kn 4 1, the molecules will almost exclusively collide with the walls of the porous structure. This case is expressed by the Knudsen diffusion number DiK given by sffiffiffiffiffiffiffiffiffiffiffi daver: 8kB T ð4Þ DiK ¼ 3 π mi with the molecular mass mi. Additionally, due to the hindered molecule movement, the bulk diffusivity in a porous medium decrease, expressed, e.g. in Ho and Webb (2006), by Dij ¼ ε=τD0ij with the tortuosity τ Z 1. To take all these wall collisions into account, the diffusion coefficient is expanded using the Bosanquet formula (Veldsink et al., 1995): 1 1 1 þ ð5Þ Deff ¼ DiK Dij The effective diffusion coefficient Deff applied to Eq. (3) leads to an extended Fickian model, which is in the following used for validations of interfacial fluxes. 2.1.2. Dusty gas model The dusty gas model (DGM) describes the porous media as a stationary phase of uniformly distributed dust particles. The flux expression in the DGM for a species i in one dimension is (Veldsink et al. (1995) and Ho and Webb (2006)):
2.1. Gas diffusion models
xi N j xj N i N i 1 ∂ðxi pÞ : ¼ RT ∂z D D ij iK jai
∑
Besides the conservation of mass for each participating component i, the general relation to describe mass transport with respect to the concentrations ci of a component i, includes convection, diffusion and sinks/sources due to chemical reactions (Bird et al., 2002): ∂ci ¼ ð∇ci vÞ ð∇Ji Þ þ Ri ∂t
ð1Þ
with the velocity v, the flux Ji and the rate of production of moles Ri of species i per unit volume. In porous media such as catalytic or gas sensor layers, only a small volume on the top of the porous layer experiences convection (described using the Brinkman equation Nield and Bejan, 2006). While only the boundary region at the top of the porous layer is affected by the convective flow, diffusion is the main transport process within the middle and the bottom region of the layer. In most cases the convective part can be neglected without loss of generality. In this work chemical reactions are not considered. With these assumptions the mass transport equation (1) simplifies to
ε ∂xi p ¼ ∇N i RT ∂t
ð2Þ
where the concentration ci is replaced by the mole fraction xi using the ideal gas law relation (ci ¼ ni =V ¼ pxi =ðRTÞ), with the porosity ε and the molar flux Ni with respect to stationary axes. The flux expresses the rate of mass transport and can be calculated by various methods, e.g. using the 1st Fickian law: Ni ¼
D0ij p RT
ð3Þ
∇xi
with the bulk molecular diffusion
D0ij .
2.1.1. Extended Fickian model Within porous media, the diffusion of gas molecules differs from bulk diffusion, expressed by D0ij for molecules of type i within molecules of type j, if the mean free path of a molecule is in the
j ¼ 1;
ð6Þ
In this equation, pressure gradients can be omitted, which corresponds to thin porous layers. For a two component system, after rearranging Eq. (6) and use of Grahams relationpbetween ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi fluxes and mole masses of the components, N 2 =N1 ¼ M 1 =M 2 , N1 can be expressed by p D12 D1K ∂x1 ð7Þ RT D12 þ D1K ax1 ∂z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with a ¼ 1 M 1 =M 2 . The insertion of Eq. (7) into Eq. (2) leads to partial differential equation (PDE) of 2nd kind and can be solved either by a PDE solver (e.g. pdepe-solver of Matlab) or by an ordinary finite difference scheme. Average pore diameter and ε required for DGM were taken from an analysis of the generated porous layers pffiffiffi while the tortuosity was calculated by τ ¼ ð 3 εÞ 1 (Millington, 1959). The DGM is taken as reference for gas diffusion calculation in a porous layer. N1 ¼
2.2. Direct simulation Monte Carlo method For very small Knudsen numbers (Kn o 0:01) the Navier–Stokes continuum model can be used to describe fluid flows. Otherwise, particle based models based on the Boltzmann equation are used. However, there is neither an analytical solution to the Boltzmann equation, nor a simple numerical solution for a partial differential equation with an integral collision term (Bird, 1994). The DSMC method circumvents these difficulties through directly simulating molecule collisions with the help of methods provided for the solution of the Boltzmann equation (e.g. the Chapman–Enskog theory or the Bhatnagar–Gross–Krook (BGK) model) to link the models to the well accepted Boltzmann equation. The essential feature of DSMC is the decoupling of the molecular motion and the intermolecular collisions over a small time step. The molecules are represented as real physical
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
71
quantities, however, the number of molecules in the calculation can be reduced through so called DSMC-particles that represent a fixed number of gas molecules. The number of real gas molecules with respect to DSMC-particles is the scaling factor f. The DSMC-particles undergo collisions according to a particular collision model. For instance, the hard-sphere (HS) model describes the collision of two particles as a function of the collision parameter b (Bird, 1994): χ b ¼ d12 cos ð8Þ 2 with the mean diameter of the colliding molecules d12 ¼ ðd1 þ d2 Þ=2 and the deflection angle χ. However, the scattering law given by the HS model is physically not realistic. A more realistic model is the variable-hard-sphere (VHS) model which considers an effective scattering cross section that decreases with increasing translational energy. This relation is expressed in a variable hard sphere molecule diameter dVHS (Koura and Matsumoto, 1991): rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c eω dVHS ¼ ð9Þ
π
with a cross section constant c ¼ 3π A2 ðνÞðκ =2Þω and relative collision energy e. Values for the collision integral A2 ðνÞ, the inverse power law force exponent ν, the force constant κ and the energy exponent ω can be found in Chapman and Cowling (1953) or Koura and Matsumoto (1991). The VHS model leads – in contrast to the HS model – to a physically realistic temperature dependent viscosity coefficient μ. However, the diffusion coefficient is not described accurately with this model (see Section 4). A further modification of the impact parameter is the energy dependent exponent of the cosine, α, of the deflection angle which enables a correct reproduction of the diffusion coefficient (Koura and Matsumoto (1991)): χ b ¼ dVSS cos α : ð10Þ 2 Due to the discretization of DSMC-particle collisions, simulation truncation errors can occur if the iteration step size Δt, mesh cell size a, or f are chosen incorrectly. It should be noted that only one collision per DSMC-particle and Δt is calculated and recollisions are neglected (Garcia and Wagner, 2000). Furthermore, binary collisions of DSMC-particles are only considered within each mesh cell and too low numbers of DSMC-particles per cell result in wrong quantities of collisions (Sun et al., 2011). Thus, a compromise between minimizing the simulation errors (large a; small f and Δt), minimizing computation resources (large f and Δt), and maximize the spatial resolution of the simulation (small a) must be found. The errors occurring from chosen Δt, a, and pf ffiffiffican be estimated by calculating the mean free path λ ¼ ð 2π nd2m Þ 1 (with the number density n and molecule diapffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi meter dm ) and the most probable thermal speed ν ¼ 2kT=m. A Δt dependent truncation error was shown to be negligible if Δt r λ ð8νÞ 1 (Garcia and Wagner, 2000; Hadjiconstantinou, 2000). For a and f no noticeable error can be expected as long as the a=λ ratio is kept below 0.5 (Sun et al., 2011; Alexander et al., 1998) and the number of DSMC-particles per cell is above 5 (Sun et al., 2011).
3. Simulation method In gas sensors obtained from Flame Spray Pyrolysis (FSP) (Großmann et al., 2011; Mädler et al., 2006; Sahm et al., 2004) two interdigitated electrodes are covered with a porous layer of nano-particles, as illustrated in Fig. 1. The probe gas diffuses into the porous layer where it is oxidized/reduced. Due to the chemical
Fig. 1. Scheme of a gas sensor as synthesized by FSP. The probe gas above the nanoparticle layer diffuses into the porous structure and, in case of CO, is oxidized to CO2 (see ☆). Due to the chemical reaction the resistance of the electrode bridging porous layer changes, giving rise to the sensor signal.
reaction (rate of reaction proportional to probe gas concentration) the electrical resistance of the porous layer changes. The change in resistance is measured by the two electrodes and a signal proportional to the probe gas concentration is obtained. One key aspect of sensor response is how much of the probe gas reacts and how close to the electrodes the reactions take place (Mädler et al., 2006). Examining such process on the nano-meter scale is the aim of this work. After validating the DSMC diffusion simulations inside an empty cuboid and with symmetric gas molecules, the simulation will be extended to CO in N2 diffusion and mass transport in porous nano-structures. As the thickness of the nano-particle layer is usually much smaller than the active sensor area, the analytical diffusion solutions can be approximated by a one dimensional approach. 3.1. Analytical solution The diffusion equation in one spatial variable, and with the bulk diffusion coefficient D ¼ D0ij independent of the concentration c, is given by the initial-boundary problem to be solved for the concentration function c: 8 ∂ c ¼ D∂x2 c; x A 0; l½ ; t 4 0; > > > t > < cð0; tÞ ¼ c ; t Z 0; 1 ð11Þ > ∂x cðl; tÞ ¼ 0; t Z 0; > > > : cðx; 0Þ ¼ 0; x A 0; l½: Here, l is the distance between the inlet and the end of the calculation domain, and c1 is the constant concentration held at the inlet. The system (11) can be solved using the separation of variables, i.e. first assuming the concentration c is the product cðx; tÞ ¼ XðxÞTðtÞ of functions X and T that respectively depend on the space and the time variable only, and then enforcing the initial and boundary conditions in (11). This separated-variables solution is cD ðx; tÞ ¼ c1
4c1
π
expð Dtðð2p þ 1Þπ =2lÞ2 Þ ð2p þ1Þπ x sin ; 2p þ 1 2l p¼0 1
∑
x A 0; l½; t 40:
ð12Þ
The subscript D indicates the value of the diffusion coefficient associated with the concentration function. In order to find a value of the diffusion coefficient D that results in a good match between (12) and the numerically simulated solution cnum at a fixed time t0, a least-squares optimization and minimization of the L2 distance is used as distðDÞ ¼ ∑ jcD ðxj ; t 0 Þ cnum ðxj ; t 0 Þj2 ; xj A G
D A I;
ð13Þ
72
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
where G 0; l½ is the discrete set of positions x at which the simulated solution cnum is available, and I R is the considered range of diffusion coefficients. In the numerical evaluation of the distance distðDÞ, the series in (12) is appropriately truncated. It is expected that dist(D) attains a global minimum at a value of the diffusion coefficient D that results in a good representation of the simulation results. 3.2. DSMC simulation The Direct Simulation Monte Carlo simulations are executed with the open source software OpenFoam version 2.01 and the implemented solver dsmcFoam. This package was extended with the VSS binary collision model using Eq. (10) and the VSS collision cross section sVSS based on the viscosity coefficient μ (Bird, 1994): sVSS ¼ π
5ðα þ 1Þðα þ 2Þðm=π Þ0:5 ðkB T ref Þω
ð14Þ
16αΓ ð9=2 ωÞμEω 0:5
with the translational kinetic energy E, molecular mass m, and the gamma function Γ . Simulations were conducted inside a cuboid with height l and cubic cross section with length z. A reflecting wall with Maxwellian thermal boundary condition was implemented at height variable x¼l while an open boundary was set at height x¼0. All other boundaries were cyclic for a 2D diffusion approach. Pristine gas was initialized within the volume at t¼0 while a desired gas concentration c1 was set at the inlet for t 40. A constant gas molecule number density of 2.68666 1025 m 3 was set within the volume as well as the open boundary (inlet) representing the standard number density at 273.15 K and 1 atm (Bird, 1994). Values of a¼ 8 nm, Δt ¼ 5 10 12 s, f¼1, and N2/CO parameters shown in Table 1 were utilized for all simulations. Larger a would limit the structure resolution of the porous layer while larger f would result in too low DSMC-particle per cell values. Thus, a direct numerical representation of all molecules in the domain has been done with respect to the scaling factor f¼1. With the chosen parameters the number of DSMC-particles per patch and gas type that are introduced each time step (particle flux accumulator; pFA) is much smaller than 1. This causes a pulsed molecule inflow with the implemented dsmcFoam inflow model (see Fig. S1 in the Supporting Information for details) as the residual of the integer number of generated DSMC-particles is added to pFA of the following iteration step. If pFA is much smaller than 1, too large or negative pFA values are calculated by this method. Therefore, the code was adjusted so that pFA is not influenced by the previous residual. Regarding possible truncation errors the chosen parameters result in Δt λ ð30νÞ 1 , a λ 1 o 7:5 1 , and an average DSMC-particle number per cell 4 13. Hence, truncation errors can be expected to be neglectable small in all cases. The spatial gas composition within the cuboid was estimated from z2 10 nm volumes assuming ideal gas. All simulations were performed three times and averaged. To generate porous layers as formed during FSP gas sensor synthesis, the raw data from Mädler et al. (2006) were utilized. Single mesh cells at the positions of the layer's nano-particles (NP) were removed from the empty cuboid with OpenFoam's tool snappyHexMesh (sHM). At the gas-NP interface a Maxwellian thermal boundary condition was implemented. Typically,
snappyHexMesh (sHM) refines the removed cubes towards spheres by decreasing a and deforming the rectangular mesh cells. Even though this procedure would result in more realistic NP shapes, a DSMC-particle number below 5 should be avoided (Sun et al., 2011) to prevent inaccurate binary collision numbers. Therefore, the refinement of sHM was disabled resulting in cubic NPs with length a. The validity of this structure simplification was verified by applying an image based method developed by Yang et al. (2009) to calculate pore size distributions. Binary images of slices through the porous layer perpendicular to x were generated every 50 nm and with a resolution of 1120 1120 pixel (0.25 nm per pixel). Calculated pore sizes were averaged over all slices. The layer's porosity ε was calculated by dividing the volume removed by sHM with the total volume of the cuboid.
4. Results 4.1. Evaluation of VHS vs. VSS Self-diffusion (ε ¼100%) of N2 in N2 was computed with the VHS and VSS model and compared to the analytical solution (Eq. (12)). A cuboid with l ¼4000 nm and square cross section z¼100 nm was chosen resulting in a Kn¼ 0.017. Due to the large inlet-wall distance and investigated case of self-diffusion (i.e. wallgas phase equilibrium with N2 expected before simulation starts) effects such gas adsorption can be neglected. Fig. 2 shows the N2
Fig. 2. Self-diffusion of 20 vol.% N2 utilizing two different binary collision models for the DSMC simulation. Four different time steps are shown and compared to the analytical solution of Fick's law of diffusion with a diffusion coefficient of 1.783 10 5 m2 s 1.
Table 1 Summary of parameters used for the DSMC simulations. All values are given for 273 K and 1 atm (Bird, 1994; Chapman and Cowling, 1953.) Parameter
Molecular mass
Molecular diameter
Viscosity index
Viscosity coefficient
Scattering factor
Symbol Unit Collision model
m kg VHS/VSS
d m VHS
ω – VHS/VSS
μ N m s1 VSS
α – VSS
N2 CO
4.65 10 26 4.65 10 26
4.17 10 10 –
0.74 0.73
1.656 10 5 1.635 10 5
1.36 1.49
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
isotope concentration calculated with the three methods at three different time steps. The VHS model results in too low N2 concentrations even after 40 ns when compared to the analytical results. The difference further increases with diffusion time indicating that the VHS model results in a too low diffusion coefficient. The concentrations computed with the VSS model match well the analytical solution. As numerous literature constants such as μ, α, ω, and D were used for the simulations as well as for the analytical solution, the obtained agreement is remarkably good. A quantitative comparison between the VHS and VSS model by means of least-squares optimization (Eq. (13)) leads to D ¼1.33 10 5 m2 s 1 (not shown) and D ¼1.62 10 5 m2s 1 (Fig. 4 inset), respectively. The reported D for N2 self-diffusion is 1.783 10 5 m2 s 1 (Massman, 1998) and shows the significant improvement of DSMC diffusion simulations using the VSS model. The validity of the DSMC method and VSS model was further evaluated by comparison with the well established semi-infinite diffusion equation (Cussler, 1997):
73
into the x 4 4000 nm region. Consequently, the gas concentration close to 4000 nm is lower in the latter case. Again, the DSMC solver combined with the VSS model results in good agreement to the analytic solution for cinf: .
4.2. Evaluation of Kn influence
with the inlet gas concentration c1 and error function erf. As OpenFOAM's DSMC solver contains only a single inflow boundary model (in- and outflow not possible) and the implementation of an additional boundary model exceeds the scope of this work, cinf: was compared to an empty cuboid with l ¼10,000 nm. Only results are shown where all regarded DSMC-particles are more than 1000 nm away from the reflecting wall at x ¼l. For clarity only the x ¼2500–4000 nm region is shown and results are compared to DSMC results of diffusion in a 4000 nm cuboid (i.e. wall at x ¼4000 nm). As can be seen in Fig. 3 the simulation and analytic results are identical for small t due to non-existences of gas molecule-wall interactions. With increasing t the DSMC-particles start to be reflected from the wall in the 4000 nm cuboid. Due to the absence of a wall in the 10,000 nm cuboid, the gas can diffuse
The influence of the Knudsen number on the diffusion is shown in Fig. 4. With increasing Kn the N2 self-diffusion is slowed down indicating a continuum to slip flow transition. The influence of particle-wall interactions increases with decreasing wall-inlet distance and thus larger Kn number. In case of large wall-inlet distances the reflecting wall has negligible influence on the diffusion and the simulation results converge towards the analytical solution for continuum regime. The low inlet concentration of 18.5 vol.% in case of Kn ¼0.068 can be explained by insufficient time to increase the N2 isotope concentration within the first 10 nm to 20 vol.% (see Fig. S2 in Supporting Information for details). The inset of Fig. 4 shows the Kn number dependance on the D. Large Kn numbers result in smaller D (i.e. D¼ 1.21 10 5 m2 s 1 for Kn ¼0.068) than the one reported for the continuum regime (1.783 10 5 m2 s 1 Massman, 1998). With decreasing Kn number D increases and approaches its final value of 1.66 10 5 m2 s 1, which is slightly lower than the D used in the analytical solution (1.783 10 5 m2 s 1 Massman, 1998). One common gas of interest in sensor or catalytic applications is CO in air. Therefore, N2 was chosen as a first approximation of air and the diffusion of 1 vol.% CO was simulated. An empty cuboid with length l ¼1000 nm and square cross section with z ¼280 nm was used to evaluate the suitability of dsmcFoam for binary gas mixtures and asymmetric gas molecules such as CO. The simulation results are again qualitatively in good agreement with the analytical solution (Fig. 5). A slightly lower DSMC diffusion in case of Kn ¼0.068 can be expected due to a continuum to slip-flow transition (fitted D ¼1.51 10 5 m2 s 1 to reported D ¼1.804
Fig. 3. Self-diffusion of 20 vol.% N2 in a cuboid with l of 4000 nm and 10,000 nm and comparison with Eqs. (12) and (15), respectively. Both cases are practically identical until the DSMC-particles start to interact with the reflecting wall in the 4000 nm cuboid. As the wall at 4000 nm is absent in the 10,000 nm cuboid, the gas diffuses in the x 4 4000 nm region results in lower concentrations close to x ¼4000 nm.
Fig. 4. N2 self-diffusion within an empty cuboid with different heights l plotted 2 1 after t l ¼ 7500 s m 2 over the relative distance x~ ¼ x l . The DSMC results are compared with the analytical solution of Ficks law of diffusion and D¼ 1.783 10 5 m2 s 1. Inset: Fitted diffusion coefficient for DSMC results taking the minimum of the sum of squared residuals between the simulation and analytical solution with varying D.
x cinf: ¼ c1 c1 erf pffiffiffiffiffiffiffiffi 4Dt
ð15Þ
74
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
10 5 m2 s 1 Massman, 1998) as discussed previously for the case of N2 self-diffusion in Fig. 4.
4.3. Diffusion in porous layers Porous layers were generated with two different structure precisions, i.e. the original one with spherical NPs with 10 nm diameters and a layer simplified with 8 nm cubes (Fig. 8). Cubes with 8 nm length were chosen to gain a similar particle volume compared to the 10 nm spheres. The porosity was constant for both layers (spherical ε ¼93% and cubic ε ¼92.6%). An example of the image based pore size calculations is shown in Fig. 6 (particles in black). A similar pattern in particle distribution and pore sizes is visible. Calculated pore diameters dpore reach in both cases from 8 to 110 nm. To further quantify the difference between the original and simplified layers the pore area A at each dpore , the normalized differential pore area ΔA ðlog ½Δdpore Σ AÞ 1 , and median pore diameter based on the Q2 pore size distribution
Fig. 5. Diffusion of 1 vol.% CO in N2 within an empty cuboid (l ¼ 1000 nm) plotted at different t. The simulation results are slightly slower than the analytical solution calculated with D ¼1.804 10 5 m2 s 1 as can be expected for Kn¼ 0.068 (Fig. 4).
d50;2 was calculated (Fig. 7). The normalized differential pore area shows that both layers (original and simplified) exhibit a broad peak between 13–100 nm. The sharp peak at 8 nm in the simplified structure can be explained by the structured gaps between the cubes (see Fig. S3 in Supporting Information for details). As the pore area difference between both layers around 8 nm (see Fig. 7 inset) is negligibly small, no significant influence is expected. The ratio of overall pore area to image area results in 92.8% and 92.6% for the original and simplified structure, respectively, which agrees well with the volume based porosity (93% and 92.6%). Also d50;2 of both structures are very similar (58 nm for original and 60 nm for simplified; see inset Fig. 7) and agree with experimental pore size distributions of FSP layers (d50;3 ¼ 70 nm for Q3 pore size distribution of as-synthesized layer Schopf et al., 2013). As the wallmolecule collisions in dsmcFoam are exclusively diffusive (i.e. reflection angle independent from impact angle) and ε as well as
Fig. 7. Normalized differential pore area of the original and simplified porous layer. Overall pore to image area ratios of 92.8% and 92.6% for the original and simplified structures can be calculated as well as d50,2 values of 58 nm and 60 nm (see inset). The peak at 8 nm in the normalized differential pore area for the simplified structure can be explained by the gaps between the 8 nm cubes (see Fig. S3 in Supporting Information for details).
Fig. 6. Color-scaled pore sizes within a slice through the original (a) and simplified (b) layers. Both structures show good agreement in pore sizes and particle (shown in black) distribution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
75
Fig. 8. Example for porous structure with 8 nm cubes taken for DSMC diffusion simulations. Shown is the FSP layer with a porosity of 92.6% and d50;2 of 60 nm.
Fig. 10. CO concentration profile generated with DSCM simulations in layers with two different porosities and comparison to a isotropic layer as well as DGM calculations. Results are shown after 5 ns and 20 ns.
Fig. 9. Normalized differential pore area of two layers as formed by FSP (ε of 92.6% and 98.1%) and comparison with a isotropic structure (ε¼ 92.6%). The overall pore to image area ratios result in 92.6% (FSP ε ¼92.6%), 92.6% (isotropic ε¼ 92.6%), and 98% (FSP ε ¼98.1%). Calculated d50;2 values are 60 nm, 41 nm, and 144 nm, respectively (see inset).
d50;2 of the original and simplified layers showed good agreement, the structure simplified with 8 nm cubes was taken for the following diffusion simulations (Fig. 8). Two additional porous structures where generated for a more comprehensive inside into layer parameters influencing the diffusion of gases: (1) An isotropic layer (randomly distributed particles) with same particle size and porosity as above and (2) a porous layer as formed during FSP but with higher porosity (98.1%). The former one is unrealistic from a mechanical point of view (floating nano-particles) but can be utilized for better evaluation and comparison with DGM. The isotropic layer exhibits a volume and area based porosity of 92.6%, which matches the FSP layer with 92.6%. The pore distribution on the other hand is shifted towards lower values (Fig. 9), which is also reflected by the smaller d50;2 (41 nm compared to 60 nm; inset Fig. 9). The decrease in d50;2 can be expected as the larger pores caused by the aggregate structure of the nano-particles (Fig. 6 and S3 in Supporting Information; see Ref. Mädler et al., 2006 for more details) are absent in the isotropic layer. An increase in porosity from 92.6% to 98.1% results in lager pores and an d50;2 of 144 nm (Fig. 9). The diffusion through the three layers characterized in Fig. 9 was simulated using DSMC and results are compared to the DGM model (Fig. 10). Note that CO adsorption was not considered in the
present study. Extension from not CO adsorbing layers to e.g. SnO2 supported Pt may result in deviations. In this case adsorption should be considered in the wall interaction scheme of the model. The DSMC results for all three films agree well to the results generated with the DGM. A comparison between the two FSP layers shows that a smaller porosity result in slower diffusion. This can be expected as a larger volume of nano-particles in the layer increases the number of particle-molecule collisions and hereby hinders its diffusion. This effect is even more pronounced with respect to the empty cuboid (i.e. ε ¼100%; Fig. 5). At t¼5 ns the CO concentration between x ¼800–1000 nm is close to 0 vol.% independent of the porosity (Fig. 10) while it already reached about 0.03 vol.% in case of the empty volume (Fig. 5). Comparing the FSP and isotropic layer with ε ¼92.6% shows the influence of pore size on the CO diffusion. Larger pores result in a larger area for the gas molecules to go through. However, the average pore size and tortuosity of inhomogeneous layers (i.e. FSP layers) can not be derived simply from the porosity. The same tortuosity was utilized for the DGM calculations of the FSP and isotropic layer as the porosity is identical (i.e. τ ¼ ðε1=3 Þ 1 Millington, 1959). It can be expected from the different pore size distributions (Fig. 9) that this assumption is inaccurate in case of the FSP layer, which explains the better match of the diffusion in the isotropic layer in Fig. 10. Thus DGM forfeits its accuracy with increasing discrepancies from the here considered isotropic layer with physically inadequate floating particles. Also, DSMC solvers require no pore diameter, tortuosity, or porosity information of the porous structure in contrast to other theoretical methods such as DGM. In order to gain such structural information time intensive and laborious structure analysis methods are required, e.g. the image based pore size calculations applied in this study. Even though DGM requires much less computational time DSMC is more efficient in terms of computational resources when the structure analysis process is considered in the overall workload. Thus, DSMC can be easily applied to simulate diffusion in graded or inhomogeneous structures where other theoretical methods reach certain limitations. With the correct mass transport verified e.g. chemical reactions could be implemented in DSMC (Bird, 1994). Investigations such as
76
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
ε, tortuosity, particle size, surface area, or particle layer thickness on sensor response time, gas selectivity or catalytic activity would be of significant scientific relevance. 5. Conclusions An open-source DSMC code was extended by the VSS binary collision model to simulate gas diffusion in nano-scaled porous layers. A pulsed DSMC-particle initialization in case of a particle flux accumulator smaller than 1 was resolved by modifying the existing inlet model. The simulation results were verified with the analytical diffusion equation and DGM. It was shown that the implementation of the VSS model significantly enhances the accuracy of DSMC diffusion simulations and a pulsed inlet condition can be avoided. Good agreement was observed between the DGM and DSMC simulations in highly porous isotropic layers. However, no prior knowledge about Kn numbers, porosity, or tortuosity as for the analytical equation or DGM are required for DSMC. Furthermore, DSMC is not limited to homogenous structures. With accurate mass transport on the nano-scale by DSMC simulations, further extensions such as gas adsorption and chemical reactions are possible. Acknowledgments LM would like to thank the German Research Foundation (DFG) for funding this project under Grant MA 3333/2-1. JD is thankful to Sven Schopf for helpful discussions regarding pore size distributions. Appendix A. Supplementary data Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ces.2013.10.038. References Alexander, F.J., Garcia, A.L., Alder, B.J., 1998. Cell size dependence of transport coefficients in stochastic particle algorithms. Phys. Fluids 10, 1540–1542.
Bird, G.A., 1994. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Clarendon Press; Oxford. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 2002. Transport Phenomena. Wiley. Chapman, S., Cowling, T.G., 1953. The Mathematical Theory of Non-Uniform Gases. Cambridge University Press, London. Cunningham, R.E., Williams, R.J.J., 1980. Diffusion in Gases and Porous Media. Plenum Press, New York. Cussler, E.L., 1997. Diffusion: Mass Transfer in Fluid Systems. Cambridge University Press. Garcia, A.L., Wagner, W., 2000. Time step truncation error in direct simulation Monte Carlo. Phys. Fluids 12, 2621–2633. Großmann, K., Kovács, K.E., Pham, D.K., Mädler, L., Barsan, N., Weimar, U., 2011. Enhancing performance of FSP SnO2-based gas sensors through Sb-doping and Pd-functionalization. Sens. Actuat. B: Chem. 158, 388–392. Hadjiconstantinou, N.G., 2000. Analysis of discretization in the direct simulation Monte Carlo. Phys. Fluids 12, 2634–2638. Ho, C.K., Webb, S.W., 2006. Gas Transport in Porous Media. Springer, Dordrecht. Koura, K., Matsumoto, H., 1991. Variable soft sphere molecular model for inversepower-law or Lennard–Jones potential. Phys. Fluids A 3, 2459–2465. Mädler, L., Roessler, A., Pratsinis, S.E., Sahm, T., Gurlo, A., Barsan, N., Weimar, U., 2006. Direct formation of highly porous gas-sensing films by in situ thermophoretic deposition of flame-made Pt/SnO2 nanoparticles. Sens. Actuat. B: Chem. 114, 283–295. Mädler, L., Lall, A.A., Friedlander, S.K., 2006. One-step aerosol synthesis of nanoparticle agglomerate films: simulation of film porosity and thickness. Nanotechnology 17, 4783–4795. Massman, W.J., 1998. A review of the molecular diffusivities of H2O, CO2, CH4, CO, O3, SO2, NH3, N2O, NO, and NO2 in air, O2 and N2 near STP. Atmos. Environ. 32, 1111–1127. Millington, R.J., 1959. Gas diffusion in porous media. Science 130, 100–102. Nield, D.A., Bejan, A., 2006. Convection in Porous Media. Springer, New York. Petropoulos, J.H., Papadokostaki, K.G., 2012. May the Knudsen equation be legitimately, or at least usefully, applied to dilute adsorbable gas flow in mesoporous media? Chem. Eng. Sci. 68, 392–400. Sahm, T., Mädler, L., Gurlo, A., Barsan, N., Pratsinis, S.E., Weimar, U., 2004. Flame spray synthesis of tin dioxide nanoparticles for gas sensing. Sens. Actuat. B: Chem. 98, 148–153. Schopf, S.O., Salameh, S., Mädler, L., 2013. Transfer of highly porous nanoparticle layers to various substrates through mechanical compression. Nanoscale 5, 3764–3772. Succi, S., 2001. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Clarendon Press, Oxford. Sun, Z.-X., Tang, Z., He, Y.-L., Tao, W.-Q., 2011. Proper cell dimension and number of particles per cell for DSMC. Comput. Fluids 50, 1–9. Veldsink, J.W., van Damme, R.M.J., Versteeg, G.F., van Swaaij, W.P.M., 1995. The use of the dusty-gas model for the description of mass transport with chemical reaction in porous media. Chem. Eng. J. 57, 115–125. Yang, Z., Peng, X.-F., Lee, D.-J., Chen, M.-Y., 2009. An image-based method for obtaining pore-size distribution of porous media. Environ. Sci. Technol. 43, 3248–3253.
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: www.elsevier.com/locate/ces
Simulation of gas diffusion in highly porous nanostructures by direct simulation Monte Carlo Jochen A.H. Dreyer a,b, Norbert Riefler b, Georg R. Pesch b, Mirza Karamehmedović b, Udo Fritsching b, Wey Yang Teoh a, Lutz Mädler b,n a b
Clean Energy and Nanotechnology (CLEAN) Laboratory, School of Energy and Environment, City University of Hong Kong, Hong Kong SAR Foundation Institute of Materials Science (IWT), Department of Production Engineering, University of Bremen, Germany
H I G H L I G H T S
G R A P H I C A L
A B S T R A C T
Direct Simulation Monte Carlo of gas diffusion with the open source solver dsmcFOAM. Improved accuracy by implementation of the variable soft sphere model. Correction of the particle flux accumulator to avoid pulsed gas initialization. CO diffusion in porous anisotropic nanostructures as formed with FSP. Base for describing adsorption, desorption and chemical reactions.
art ic l e i nf o
a b s t r a c t
Article history: Received 4 September 2013 Received in revised form 17 October 2013 Accepted 22 October 2013 Available online 6 November 2013
A Direct Simulation Monte Carlo (DSMC) method is utilized to simulate gas diffusion in nanoscaled highly porous layers. An open source solver has been extended with the variable soft sphere (VSS) binary collision model and the inflow boundary model was adjusted for small numbers of DSMC particle initialization. Comparison with the analytical diffusion equation illustrate the improvement of the VSS model compared to the variable hard sphere model (VHS). Subsequently, several highly porous particle layers (gas sensors synthesized by flame spray pyrolysis and isotropic layers) build up by 10 nm particles have been investigated. Results for DSMC gas diffusion in the porous structures are in agreement with the well established dusty gas model (DGM). However, while DGM requires measurements or estimations of pore sizes, porosity, and tortuosity and furthermore is limited to homogenous layers, the present contribution shows significant advantages of DSMC in describing gas diffusion in nonisotropic porous structures. & 2013 Elsevier Ltd. All rights reserved.
Keywords: Nanoparticle aggregate films Gas diffusion DSMC Dusty gas model High porosity Gas sensor
1. Introduction Diffusion of gas molecules into porous structures is a phenomenon which occurs in many environmental or technical systems such as gas sensors and catalyst beds. The understanding of macroscopic effects in such systems often requires modeling at the mesoscopic scale. To describe rarefied gases or systems with
n
Corresponding author. Tel.: þ 49 421 218 51200; fax: þ 49 421 218 51211. E-mail address: [email protected] (L. Mädler).
0009-2509/$ - see front matter & 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.ces.2013.10.038
small pores (Knudsen number Kn 4 0:1), e.g. mesoporous systems (pore-diameter between 2…50 nm), kinetic theory has to be applied through solving the Boltzmann equation. Early approximations of this equation for pure diffusion without walls led (besides the Boltzmann H-theorem) to the Chapman–Enskog perturbation method, followed by the Grad–Zhdanov moment technique. The linearization of the Boltzmann equation was introduced later by Bhatnagar, Gross and Krook (Cunningham and Williams (1980)) and serves as the basis of many lattice Boltzmann model (LBM) simulation methods (Succi, 2001). In contrast, the dusty gas model approach (DGM) uses the
70
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
non-linear Chapman–Enskog or the Grad–Zhdanov constitutive equations of diffusion, where wall effects are considered by collisions of the gas molecules with stationary rigid spheres (Cunningham and Williams, 1980). Another approach to calculate diffusive and other kinetic properties is the direct tracking, i.e. calculation of tracks and collisions of gas molecules. Molecular dynamics (MD) is an often applied method where each atom or molecule is individually tracked, however, with the expense of very high computational costs if larger problems are considered. A good compromise between the required computing power and problem size is the Direct Simulation Monte Carlo method (DSMC) (Bird, 1994). In this method, the molecules of a volume cell are united into different parcels. The translation and the collision process are decoupled in every time step. In the collision process, only a part of the molecules are chosen according to their collision cross section. These collision partners represent statistically the ensemble of all molecules within the simulation volume. This work describes the implementation of the variable-softsphere (VSS) collision model (Koura and Matsumoto, 1991) in an existing DSMC solver of OpenFOAM. The results are compared with analytical solutions of the diffusion equation. Beyond that, the diffusion processes into gas sensor layers are investigated based on this implementation and compared to the well established DGM method.
2. Theoretical models
size range of the pore diameter. This case occurs for Knudsen numbers Kn Z 0:01 with Kn ¼ λ=daver: , where daver: is the average pore diameter and λ the mean free path. In case of non-adsorbing gas molecules (Petropoulos and Papadokostaki, 2012) and for Kn 4 1, the molecules will almost exclusively collide with the walls of the porous structure. This case is expressed by the Knudsen diffusion number DiK given by sffiffiffiffiffiffiffiffiffiffiffi daver: 8kB T ð4Þ DiK ¼ 3 π mi with the molecular mass mi. Additionally, due to the hindered molecule movement, the bulk diffusivity in a porous medium decrease, expressed, e.g. in Ho and Webb (2006), by Dij ¼ ε=τD0ij with the tortuosity τ Z 1. To take all these wall collisions into account, the diffusion coefficient is expanded using the Bosanquet formula (Veldsink et al., 1995): 1 1 1 þ ð5Þ Deff ¼ DiK Dij The effective diffusion coefficient Deff applied to Eq. (3) leads to an extended Fickian model, which is in the following used for validations of interfacial fluxes. 2.1.2. Dusty gas model The dusty gas model (DGM) describes the porous media as a stationary phase of uniformly distributed dust particles. The flux expression in the DGM for a species i in one dimension is (Veldsink et al. (1995) and Ho and Webb (2006)):
2.1. Gas diffusion models
xi N j xj N i N i 1 ∂ðxi pÞ : ¼ RT ∂z D D ij iK jai
∑
Besides the conservation of mass for each participating component i, the general relation to describe mass transport with respect to the concentrations ci of a component i, includes convection, diffusion and sinks/sources due to chemical reactions (Bird et al., 2002): ∂ci ¼ ð∇ci vÞ ð∇Ji Þ þ Ri ∂t
ð1Þ
with the velocity v, the flux Ji and the rate of production of moles Ri of species i per unit volume. In porous media such as catalytic or gas sensor layers, only a small volume on the top of the porous layer experiences convection (described using the Brinkman equation Nield and Bejan, 2006). While only the boundary region at the top of the porous layer is affected by the convective flow, diffusion is the main transport process within the middle and the bottom region of the layer. In most cases the convective part can be neglected without loss of generality. In this work chemical reactions are not considered. With these assumptions the mass transport equation (1) simplifies to
ε ∂xi p ¼ ∇N i RT ∂t
ð2Þ
where the concentration ci is replaced by the mole fraction xi using the ideal gas law relation (ci ¼ ni =V ¼ pxi =ðRTÞ), with the porosity ε and the molar flux Ni with respect to stationary axes. The flux expresses the rate of mass transport and can be calculated by various methods, e.g. using the 1st Fickian law: Ni ¼
D0ij p RT
ð3Þ
∇xi
with the bulk molecular diffusion
D0ij .
2.1.1. Extended Fickian model Within porous media, the diffusion of gas molecules differs from bulk diffusion, expressed by D0ij for molecules of type i within molecules of type j, if the mean free path of a molecule is in the
j ¼ 1;
ð6Þ
In this equation, pressure gradients can be omitted, which corresponds to thin porous layers. For a two component system, after rearranging Eq. (6) and use of Grahams relationpbetween ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi fluxes and mole masses of the components, N 2 =N1 ¼ M 1 =M 2 , N1 can be expressed by p D12 D1K ∂x1 ð7Þ RT D12 þ D1K ax1 ∂z pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi with a ¼ 1 M 1 =M 2 . The insertion of Eq. (7) into Eq. (2) leads to partial differential equation (PDE) of 2nd kind and can be solved either by a PDE solver (e.g. pdepe-solver of Matlab) or by an ordinary finite difference scheme. Average pore diameter and ε required for DGM were taken from an analysis of the generated porous layers pffiffiffi while the tortuosity was calculated by τ ¼ ð 3 εÞ 1 (Millington, 1959). The DGM is taken as reference for gas diffusion calculation in a porous layer. N1 ¼
2.2. Direct simulation Monte Carlo method For very small Knudsen numbers (Kn o 0:01) the Navier–Stokes continuum model can be used to describe fluid flows. Otherwise, particle based models based on the Boltzmann equation are used. However, there is neither an analytical solution to the Boltzmann equation, nor a simple numerical solution for a partial differential equation with an integral collision term (Bird, 1994). The DSMC method circumvents these difficulties through directly simulating molecule collisions with the help of methods provided for the solution of the Boltzmann equation (e.g. the Chapman–Enskog theory or the Bhatnagar–Gross–Krook (BGK) model) to link the models to the well accepted Boltzmann equation. The essential feature of DSMC is the decoupling of the molecular motion and the intermolecular collisions over a small time step. The molecules are represented as real physical
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
71
quantities, however, the number of molecules in the calculation can be reduced through so called DSMC-particles that represent a fixed number of gas molecules. The number of real gas molecules with respect to DSMC-particles is the scaling factor f. The DSMC-particles undergo collisions according to a particular collision model. For instance, the hard-sphere (HS) model describes the collision of two particles as a function of the collision parameter b (Bird, 1994): χ b ¼ d12 cos ð8Þ 2 with the mean diameter of the colliding molecules d12 ¼ ðd1 þ d2 Þ=2 and the deflection angle χ. However, the scattering law given by the HS model is physically not realistic. A more realistic model is the variable-hard-sphere (VHS) model which considers an effective scattering cross section that decreases with increasing translational energy. This relation is expressed in a variable hard sphere molecule diameter dVHS (Koura and Matsumoto, 1991): rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c eω dVHS ¼ ð9Þ
π
with a cross section constant c ¼ 3π A2 ðνÞðκ =2Þω and relative collision energy e. Values for the collision integral A2 ðνÞ, the inverse power law force exponent ν, the force constant κ and the energy exponent ω can be found in Chapman and Cowling (1953) or Koura and Matsumoto (1991). The VHS model leads – in contrast to the HS model – to a physically realistic temperature dependent viscosity coefficient μ. However, the diffusion coefficient is not described accurately with this model (see Section 4). A further modification of the impact parameter is the energy dependent exponent of the cosine, α, of the deflection angle which enables a correct reproduction of the diffusion coefficient (Koura and Matsumoto (1991)): χ b ¼ dVSS cos α : ð10Þ 2 Due to the discretization of DSMC-particle collisions, simulation truncation errors can occur if the iteration step size Δt, mesh cell size a, or f are chosen incorrectly. It should be noted that only one collision per DSMC-particle and Δt is calculated and recollisions are neglected (Garcia and Wagner, 2000). Furthermore, binary collisions of DSMC-particles are only considered within each mesh cell and too low numbers of DSMC-particles per cell result in wrong quantities of collisions (Sun et al., 2011). Thus, a compromise between minimizing the simulation errors (large a; small f and Δt), minimizing computation resources (large f and Δt), and maximize the spatial resolution of the simulation (small a) must be found. The errors occurring from chosen Δt, a, and pf ffiffiffican be estimated by calculating the mean free path λ ¼ ð 2π nd2m Þ 1 (with the number density n and molecule diapffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi meter dm ) and the most probable thermal speed ν ¼ 2kT=m. A Δt dependent truncation error was shown to be negligible if Δt r λ ð8νÞ 1 (Garcia and Wagner, 2000; Hadjiconstantinou, 2000). For a and f no noticeable error can be expected as long as the a=λ ratio is kept below 0.5 (Sun et al., 2011; Alexander et al., 1998) and the number of DSMC-particles per cell is above 5 (Sun et al., 2011).
3. Simulation method In gas sensors obtained from Flame Spray Pyrolysis (FSP) (Großmann et al., 2011; Mädler et al., 2006; Sahm et al., 2004) two interdigitated electrodes are covered with a porous layer of nano-particles, as illustrated in Fig. 1. The probe gas diffuses into the porous layer where it is oxidized/reduced. Due to the chemical
Fig. 1. Scheme of a gas sensor as synthesized by FSP. The probe gas above the nanoparticle layer diffuses into the porous structure and, in case of CO, is oxidized to CO2 (see ☆). Due to the chemical reaction the resistance of the electrode bridging porous layer changes, giving rise to the sensor signal.
reaction (rate of reaction proportional to probe gas concentration) the electrical resistance of the porous layer changes. The change in resistance is measured by the two electrodes and a signal proportional to the probe gas concentration is obtained. One key aspect of sensor response is how much of the probe gas reacts and how close to the electrodes the reactions take place (Mädler et al., 2006). Examining such process on the nano-meter scale is the aim of this work. After validating the DSMC diffusion simulations inside an empty cuboid and with symmetric gas molecules, the simulation will be extended to CO in N2 diffusion and mass transport in porous nano-structures. As the thickness of the nano-particle layer is usually much smaller than the active sensor area, the analytical diffusion solutions can be approximated by a one dimensional approach. 3.1. Analytical solution The diffusion equation in one spatial variable, and with the bulk diffusion coefficient D ¼ D0ij independent of the concentration c, is given by the initial-boundary problem to be solved for the concentration function c: 8 ∂ c ¼ D∂x2 c; x A 0; l½ ; t 4 0; > > > t > < cð0; tÞ ¼ c ; t Z 0; 1 ð11Þ > ∂x cðl; tÞ ¼ 0; t Z 0; > > > : cðx; 0Þ ¼ 0; x A 0; l½: Here, l is the distance between the inlet and the end of the calculation domain, and c1 is the constant concentration held at the inlet. The system (11) can be solved using the separation of variables, i.e. first assuming the concentration c is the product cðx; tÞ ¼ XðxÞTðtÞ of functions X and T that respectively depend on the space and the time variable only, and then enforcing the initial and boundary conditions in (11). This separated-variables solution is cD ðx; tÞ ¼ c1
4c1
π
expð Dtðð2p þ 1Þπ =2lÞ2 Þ ð2p þ1Þπ x sin ; 2p þ 1 2l p¼0 1
∑
x A 0; l½; t 40:
ð12Þ
The subscript D indicates the value of the diffusion coefficient associated with the concentration function. In order to find a value of the diffusion coefficient D that results in a good match between (12) and the numerically simulated solution cnum at a fixed time t0, a least-squares optimization and minimization of the L2 distance is used as distðDÞ ¼ ∑ jcD ðxj ; t 0 Þ cnum ðxj ; t 0 Þj2 ; xj A G
D A I;
ð13Þ
72
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
where G 0; l½ is the discrete set of positions x at which the simulated solution cnum is available, and I R is the considered range of diffusion coefficients. In the numerical evaluation of the distance distðDÞ, the series in (12) is appropriately truncated. It is expected that dist(D) attains a global minimum at a value of the diffusion coefficient D that results in a good representation of the simulation results. 3.2. DSMC simulation The Direct Simulation Monte Carlo simulations are executed with the open source software OpenFoam version 2.01 and the implemented solver dsmcFoam. This package was extended with the VSS binary collision model using Eq. (10) and the VSS collision cross section sVSS based on the viscosity coefficient μ (Bird, 1994): sVSS ¼ π
5ðα þ 1Þðα þ 2Þðm=π Þ0:5 ðkB T ref Þω
ð14Þ
16αΓ ð9=2 ωÞμEω 0:5
with the translational kinetic energy E, molecular mass m, and the gamma function Γ . Simulations were conducted inside a cuboid with height l and cubic cross section with length z. A reflecting wall with Maxwellian thermal boundary condition was implemented at height variable x¼l while an open boundary was set at height x¼0. All other boundaries were cyclic for a 2D diffusion approach. Pristine gas was initialized within the volume at t¼0 while a desired gas concentration c1 was set at the inlet for t 40. A constant gas molecule number density of 2.68666 1025 m 3 was set within the volume as well as the open boundary (inlet) representing the standard number density at 273.15 K and 1 atm (Bird, 1994). Values of a¼ 8 nm, Δt ¼ 5 10 12 s, f¼1, and N2/CO parameters shown in Table 1 were utilized for all simulations. Larger a would limit the structure resolution of the porous layer while larger f would result in too low DSMC-particle per cell values. Thus, a direct numerical representation of all molecules in the domain has been done with respect to the scaling factor f¼1. With the chosen parameters the number of DSMC-particles per patch and gas type that are introduced each time step (particle flux accumulator; pFA) is much smaller than 1. This causes a pulsed molecule inflow with the implemented dsmcFoam inflow model (see Fig. S1 in the Supporting Information for details) as the residual of the integer number of generated DSMC-particles is added to pFA of the following iteration step. If pFA is much smaller than 1, too large or negative pFA values are calculated by this method. Therefore, the code was adjusted so that pFA is not influenced by the previous residual. Regarding possible truncation errors the chosen parameters result in Δt λ ð30νÞ 1 , a λ 1 o 7:5 1 , and an average DSMC-particle number per cell 4 13. Hence, truncation errors can be expected to be neglectable small in all cases. The spatial gas composition within the cuboid was estimated from z2 10 nm volumes assuming ideal gas. All simulations were performed three times and averaged. To generate porous layers as formed during FSP gas sensor synthesis, the raw data from Mädler et al. (2006) were utilized. Single mesh cells at the positions of the layer's nano-particles (NP) were removed from the empty cuboid with OpenFoam's tool snappyHexMesh (sHM). At the gas-NP interface a Maxwellian thermal boundary condition was implemented. Typically,
snappyHexMesh (sHM) refines the removed cubes towards spheres by decreasing a and deforming the rectangular mesh cells. Even though this procedure would result in more realistic NP shapes, a DSMC-particle number below 5 should be avoided (Sun et al., 2011) to prevent inaccurate binary collision numbers. Therefore, the refinement of sHM was disabled resulting in cubic NPs with length a. The validity of this structure simplification was verified by applying an image based method developed by Yang et al. (2009) to calculate pore size distributions. Binary images of slices through the porous layer perpendicular to x were generated every 50 nm and with a resolution of 1120 1120 pixel (0.25 nm per pixel). Calculated pore sizes were averaged over all slices. The layer's porosity ε was calculated by dividing the volume removed by sHM with the total volume of the cuboid.
4. Results 4.1. Evaluation of VHS vs. VSS Self-diffusion (ε ¼100%) of N2 in N2 was computed with the VHS and VSS model and compared to the analytical solution (Eq. (12)). A cuboid with l ¼4000 nm and square cross section z¼100 nm was chosen resulting in a Kn¼ 0.017. Due to the large inlet-wall distance and investigated case of self-diffusion (i.e. wallgas phase equilibrium with N2 expected before simulation starts) effects such gas adsorption can be neglected. Fig. 2 shows the N2
Fig. 2. Self-diffusion of 20 vol.% N2 utilizing two different binary collision models for the DSMC simulation. Four different time steps are shown and compared to the analytical solution of Fick's law of diffusion with a diffusion coefficient of 1.783 10 5 m2 s 1.
Table 1 Summary of parameters used for the DSMC simulations. All values are given for 273 K and 1 atm (Bird, 1994; Chapman and Cowling, 1953.) Parameter
Molecular mass
Molecular diameter
Viscosity index
Viscosity coefficient
Scattering factor
Symbol Unit Collision model
m kg VHS/VSS
d m VHS
ω – VHS/VSS
μ N m s1 VSS
α – VSS
N2 CO
4.65 10 26 4.65 10 26
4.17 10 10 –
0.74 0.73
1.656 10 5 1.635 10 5
1.36 1.49
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
isotope concentration calculated with the three methods at three different time steps. The VHS model results in too low N2 concentrations even after 40 ns when compared to the analytical results. The difference further increases with diffusion time indicating that the VHS model results in a too low diffusion coefficient. The concentrations computed with the VSS model match well the analytical solution. As numerous literature constants such as μ, α, ω, and D were used for the simulations as well as for the analytical solution, the obtained agreement is remarkably good. A quantitative comparison between the VHS and VSS model by means of least-squares optimization (Eq. (13)) leads to D ¼1.33 10 5 m2 s 1 (not shown) and D ¼1.62 10 5 m2s 1 (Fig. 4 inset), respectively. The reported D for N2 self-diffusion is 1.783 10 5 m2 s 1 (Massman, 1998) and shows the significant improvement of DSMC diffusion simulations using the VSS model. The validity of the DSMC method and VSS model was further evaluated by comparison with the well established semi-infinite diffusion equation (Cussler, 1997):
73
into the x 4 4000 nm region. Consequently, the gas concentration close to 4000 nm is lower in the latter case. Again, the DSMC solver combined with the VSS model results in good agreement to the analytic solution for cinf: .
4.2. Evaluation of Kn influence
with the inlet gas concentration c1 and error function erf. As OpenFOAM's DSMC solver contains only a single inflow boundary model (in- and outflow not possible) and the implementation of an additional boundary model exceeds the scope of this work, cinf: was compared to an empty cuboid with l ¼10,000 nm. Only results are shown where all regarded DSMC-particles are more than 1000 nm away from the reflecting wall at x ¼l. For clarity only the x ¼2500–4000 nm region is shown and results are compared to DSMC results of diffusion in a 4000 nm cuboid (i.e. wall at x ¼4000 nm). As can be seen in Fig. 3 the simulation and analytic results are identical for small t due to non-existences of gas molecule-wall interactions. With increasing t the DSMC-particles start to be reflected from the wall in the 4000 nm cuboid. Due to the absence of a wall in the 10,000 nm cuboid, the gas can diffuse
The influence of the Knudsen number on the diffusion is shown in Fig. 4. With increasing Kn the N2 self-diffusion is slowed down indicating a continuum to slip flow transition. The influence of particle-wall interactions increases with decreasing wall-inlet distance and thus larger Kn number. In case of large wall-inlet distances the reflecting wall has negligible influence on the diffusion and the simulation results converge towards the analytical solution for continuum regime. The low inlet concentration of 18.5 vol.% in case of Kn ¼0.068 can be explained by insufficient time to increase the N2 isotope concentration within the first 10 nm to 20 vol.% (see Fig. S2 in Supporting Information for details). The inset of Fig. 4 shows the Kn number dependance on the D. Large Kn numbers result in smaller D (i.e. D¼ 1.21 10 5 m2 s 1 for Kn ¼0.068) than the one reported for the continuum regime (1.783 10 5 m2 s 1 Massman, 1998). With decreasing Kn number D increases and approaches its final value of 1.66 10 5 m2 s 1, which is slightly lower than the D used in the analytical solution (1.783 10 5 m2 s 1 Massman, 1998). One common gas of interest in sensor or catalytic applications is CO in air. Therefore, N2 was chosen as a first approximation of air and the diffusion of 1 vol.% CO was simulated. An empty cuboid with length l ¼1000 nm and square cross section with z ¼280 nm was used to evaluate the suitability of dsmcFoam for binary gas mixtures and asymmetric gas molecules such as CO. The simulation results are again qualitatively in good agreement with the analytical solution (Fig. 5). A slightly lower DSMC diffusion in case of Kn ¼0.068 can be expected due to a continuum to slip-flow transition (fitted D ¼1.51 10 5 m2 s 1 to reported D ¼1.804
Fig. 3. Self-diffusion of 20 vol.% N2 in a cuboid with l of 4000 nm and 10,000 nm and comparison with Eqs. (12) and (15), respectively. Both cases are practically identical until the DSMC-particles start to interact with the reflecting wall in the 4000 nm cuboid. As the wall at 4000 nm is absent in the 10,000 nm cuboid, the gas diffuses in the x 4 4000 nm region results in lower concentrations close to x ¼4000 nm.
Fig. 4. N2 self-diffusion within an empty cuboid with different heights l plotted 2 1 after t l ¼ 7500 s m 2 over the relative distance x~ ¼ x l . The DSMC results are compared with the analytical solution of Ficks law of diffusion and D¼ 1.783 10 5 m2 s 1. Inset: Fitted diffusion coefficient for DSMC results taking the minimum of the sum of squared residuals between the simulation and analytical solution with varying D.
x cinf: ¼ c1 c1 erf pffiffiffiffiffiffiffiffi 4Dt
ð15Þ
74
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
10 5 m2 s 1 Massman, 1998) as discussed previously for the case of N2 self-diffusion in Fig. 4.
4.3. Diffusion in porous layers Porous layers were generated with two different structure precisions, i.e. the original one with spherical NPs with 10 nm diameters and a layer simplified with 8 nm cubes (Fig. 8). Cubes with 8 nm length were chosen to gain a similar particle volume compared to the 10 nm spheres. The porosity was constant for both layers (spherical ε ¼93% and cubic ε ¼92.6%). An example of the image based pore size calculations is shown in Fig. 6 (particles in black). A similar pattern in particle distribution and pore sizes is visible. Calculated pore diameters dpore reach in both cases from 8 to 110 nm. To further quantify the difference between the original and simplified layers the pore area A at each dpore , the normalized differential pore area ΔA ðlog ½Δdpore Σ AÞ 1 , and median pore diameter based on the Q2 pore size distribution
Fig. 5. Diffusion of 1 vol.% CO in N2 within an empty cuboid (l ¼ 1000 nm) plotted at different t. The simulation results are slightly slower than the analytical solution calculated with D ¼1.804 10 5 m2 s 1 as can be expected for Kn¼ 0.068 (Fig. 4).
d50;2 was calculated (Fig. 7). The normalized differential pore area shows that both layers (original and simplified) exhibit a broad peak between 13–100 nm. The sharp peak at 8 nm in the simplified structure can be explained by the structured gaps between the cubes (see Fig. S3 in Supporting Information for details). As the pore area difference between both layers around 8 nm (see Fig. 7 inset) is negligibly small, no significant influence is expected. The ratio of overall pore area to image area results in 92.8% and 92.6% for the original and simplified structure, respectively, which agrees well with the volume based porosity (93% and 92.6%). Also d50;2 of both structures are very similar (58 nm for original and 60 nm for simplified; see inset Fig. 7) and agree with experimental pore size distributions of FSP layers (d50;3 ¼ 70 nm for Q3 pore size distribution of as-synthesized layer Schopf et al., 2013). As the wallmolecule collisions in dsmcFoam are exclusively diffusive (i.e. reflection angle independent from impact angle) and ε as well as
Fig. 7. Normalized differential pore area of the original and simplified porous layer. Overall pore to image area ratios of 92.8% and 92.6% for the original and simplified structures can be calculated as well as d50,2 values of 58 nm and 60 nm (see inset). The peak at 8 nm in the normalized differential pore area for the simplified structure can be explained by the gaps between the 8 nm cubes (see Fig. S3 in Supporting Information for details).
Fig. 6. Color-scaled pore sizes within a slice through the original (a) and simplified (b) layers. Both structures show good agreement in pore sizes and particle (shown in black) distribution. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this paper.)
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
75
Fig. 8. Example for porous structure with 8 nm cubes taken for DSMC diffusion simulations. Shown is the FSP layer with a porosity of 92.6% and d50;2 of 60 nm.
Fig. 10. CO concentration profile generated with DSCM simulations in layers with two different porosities and comparison to a isotropic layer as well as DGM calculations. Results are shown after 5 ns and 20 ns.
Fig. 9. Normalized differential pore area of two layers as formed by FSP (ε of 92.6% and 98.1%) and comparison with a isotropic structure (ε¼ 92.6%). The overall pore to image area ratios result in 92.6% (FSP ε ¼92.6%), 92.6% (isotropic ε¼ 92.6%), and 98% (FSP ε ¼98.1%). Calculated d50;2 values are 60 nm, 41 nm, and 144 nm, respectively (see inset).
d50;2 of the original and simplified layers showed good agreement, the structure simplified with 8 nm cubes was taken for the following diffusion simulations (Fig. 8). Two additional porous structures where generated for a more comprehensive inside into layer parameters influencing the diffusion of gases: (1) An isotropic layer (randomly distributed particles) with same particle size and porosity as above and (2) a porous layer as formed during FSP but with higher porosity (98.1%). The former one is unrealistic from a mechanical point of view (floating nano-particles) but can be utilized for better evaluation and comparison with DGM. The isotropic layer exhibits a volume and area based porosity of 92.6%, which matches the FSP layer with 92.6%. The pore distribution on the other hand is shifted towards lower values (Fig. 9), which is also reflected by the smaller d50;2 (41 nm compared to 60 nm; inset Fig. 9). The decrease in d50;2 can be expected as the larger pores caused by the aggregate structure of the nano-particles (Fig. 6 and S3 in Supporting Information; see Ref. Mädler et al., 2006 for more details) are absent in the isotropic layer. An increase in porosity from 92.6% to 98.1% results in lager pores and an d50;2 of 144 nm (Fig. 9). The diffusion through the three layers characterized in Fig. 9 was simulated using DSMC and results are compared to the DGM model (Fig. 10). Note that CO adsorption was not considered in the
present study. Extension from not CO adsorbing layers to e.g. SnO2 supported Pt may result in deviations. In this case adsorption should be considered in the wall interaction scheme of the model. The DSMC results for all three films agree well to the results generated with the DGM. A comparison between the two FSP layers shows that a smaller porosity result in slower diffusion. This can be expected as a larger volume of nano-particles in the layer increases the number of particle-molecule collisions and hereby hinders its diffusion. This effect is even more pronounced with respect to the empty cuboid (i.e. ε ¼100%; Fig. 5). At t¼5 ns the CO concentration between x ¼800–1000 nm is close to 0 vol.% independent of the porosity (Fig. 10) while it already reached about 0.03 vol.% in case of the empty volume (Fig. 5). Comparing the FSP and isotropic layer with ε ¼92.6% shows the influence of pore size on the CO diffusion. Larger pores result in a larger area for the gas molecules to go through. However, the average pore size and tortuosity of inhomogeneous layers (i.e. FSP layers) can not be derived simply from the porosity. The same tortuosity was utilized for the DGM calculations of the FSP and isotropic layer as the porosity is identical (i.e. τ ¼ ðε1=3 Þ 1 Millington, 1959). It can be expected from the different pore size distributions (Fig. 9) that this assumption is inaccurate in case of the FSP layer, which explains the better match of the diffusion in the isotropic layer in Fig. 10. Thus DGM forfeits its accuracy with increasing discrepancies from the here considered isotropic layer with physically inadequate floating particles. Also, DSMC solvers require no pore diameter, tortuosity, or porosity information of the porous structure in contrast to other theoretical methods such as DGM. In order to gain such structural information time intensive and laborious structure analysis methods are required, e.g. the image based pore size calculations applied in this study. Even though DGM requires much less computational time DSMC is more efficient in terms of computational resources when the structure analysis process is considered in the overall workload. Thus, DSMC can be easily applied to simulate diffusion in graded or inhomogeneous structures where other theoretical methods reach certain limitations. With the correct mass transport verified e.g. chemical reactions could be implemented in DSMC (Bird, 1994). Investigations such as
76
J.A.H. Dreyer et al. / Chemical Engineering Science 105 (2014) 69–76
ε, tortuosity, particle size, surface area, or particle layer thickness on sensor response time, gas selectivity or catalytic activity would be of significant scientific relevance. 5. Conclusions An open-source DSMC code was extended by the VSS binary collision model to simulate gas diffusion in nano-scaled porous layers. A pulsed DSMC-particle initialization in case of a particle flux accumulator smaller than 1 was resolved by modifying the existing inlet model. The simulation results were verified with the analytical diffusion equation and DGM. It was shown that the implementation of the VSS model significantly enhances the accuracy of DSMC diffusion simulations and a pulsed inlet condition can be avoided. Good agreement was observed between the DGM and DSMC simulations in highly porous isotropic layers. However, no prior knowledge about Kn numbers, porosity, or tortuosity as for the analytical equation or DGM are required for DSMC. Furthermore, DSMC is not limited to homogenous structures. With accurate mass transport on the nano-scale by DSMC simulations, further extensions such as gas adsorption and chemical reactions are possible. Acknowledgments LM would like to thank the German Research Foundation (DFG) for funding this project under Grant MA 3333/2-1. JD is thankful to Sven Schopf for helpful discussions regarding pore size distributions. Appendix A. Supplementary data Supplementary data associated with this article can be found in the online version at http://dx.doi.org/10.1016/j.ces.2013.10.038. References Alexander, F.J., Garcia, A.L., Alder, B.J., 1998. Cell size dependence of transport coefficients in stochastic particle algorithms. Phys. Fluids 10, 1540–1542.
Bird, G.A., 1994. Molecular Gas Dynamics and the Direct Simulation of Gas Flows. Clarendon Press; Oxford. Bird, R.B., Stewart, W.E., Lightfoot, E.N., 2002. Transport Phenomena. Wiley. Chapman, S., Cowling, T.G., 1953. The Mathematical Theory of Non-Uniform Gases. Cambridge University Press, London. Cunningham, R.E., Williams, R.J.J., 1980. Diffusion in Gases and Porous Media. Plenum Press, New York. Cussler, E.L., 1997. Diffusion: Mass Transfer in Fluid Systems. Cambridge University Press. Garcia, A.L., Wagner, W., 2000. Time step truncation error in direct simulation Monte Carlo. Phys. Fluids 12, 2621–2633. Großmann, K., Kovács, K.E., Pham, D.K., Mädler, L., Barsan, N., Weimar, U., 2011. Enhancing performance of FSP SnO2-based gas sensors through Sb-doping and Pd-functionalization. Sens. Actuat. B: Chem. 158, 388–392. Hadjiconstantinou, N.G., 2000. Analysis of discretization in the direct simulation Monte Carlo. Phys. Fluids 12, 2634–2638. Ho, C.K., Webb, S.W., 2006. Gas Transport in Porous Media. Springer, Dordrecht. Koura, K., Matsumoto, H., 1991. Variable soft sphere molecular model for inversepower-law or Lennard–Jones potential. Phys. Fluids A 3, 2459–2465. Mädler, L., Roessler, A., Pratsinis, S.E., Sahm, T., Gurlo, A., Barsan, N., Weimar, U., 2006. Direct formation of highly porous gas-sensing films by in situ thermophoretic deposition of flame-made Pt/SnO2 nanoparticles. Sens. Actuat. B: Chem. 114, 283–295. Mädler, L., Lall, A.A., Friedlander, S.K., 2006. One-step aerosol synthesis of nanoparticle agglomerate films: simulation of film porosity and thickness. Nanotechnology 17, 4783–4795. Massman, W.J., 1998. A review of the molecular diffusivities of H2O, CO2, CH4, CO, O3, SO2, NH3, N2O, NO, and NO2 in air, O2 and N2 near STP. Atmos. Environ. 32, 1111–1127. Millington, R.J., 1959. Gas diffusion in porous media. Science 130, 100–102. Nield, D.A., Bejan, A., 2006. Convection in Porous Media. Springer, New York. Petropoulos, J.H., Papadokostaki, K.G., 2012. May the Knudsen equation be legitimately, or at least usefully, applied to dilute adsorbable gas flow in mesoporous media? Chem. Eng. Sci. 68, 392–400. Sahm, T., Mädler, L., Gurlo, A., Barsan, N., Pratsinis, S.E., Weimar, U., 2004. Flame spray synthesis of tin dioxide nanoparticles for gas sensing. Sens. Actuat. B: Chem. 98, 148–153. Schopf, S.O., Salameh, S., Mädler, L., 2013. Transfer of highly porous nanoparticle layers to various substrates through mechanical compression. Nanoscale 5, 3764–3772. Succi, S., 2001. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond. Clarendon Press, Oxford. Sun, Z.-X., Tang, Z., He, Y.-L., Tao, W.-Q., 2011. Proper cell dimension and number of particles per cell for DSMC. Comput. Fluids 50, 1–9. Veldsink, J.W., van Damme, R.M.J., Versteeg, G.F., van Swaaij, W.P.M., 1995. The use of the dusty-gas model for the description of mass transport with chemical reaction in porous media. Chem. Eng. J. 57, 115–125. Yang, Z., Peng, X.-F., Lee, D.-J., Chen, M.-Y., 2009. An image-based method for obtaining pore-size distribution of porous media. Environ. Sci. Technol. 43, 3248–3253.