Simulation Study on the Reaction-Diffusion ... - ACS Publications

Article pubs.acs.org/Langmuir

Simulation Study on the Reaction-Diffusion Coupling in Simple Pore Structures Yanping Li,†,‡,§ Mingcan Zhao,‡,∥ Chengxiang Li,*,‡ and Wei Ge*,‡,§,∥ †

School of Chemical Engineering and Technology, Tianjin University, Tianjin 300072, China State Key Laboratory of Multiphase Complex Systems, Institute of Process Engineering, Chinese Academy of Sciences, Beijing 100190, China § Collaborative Innovation Center of Chemical Science and Engineering (Tianjin), Tianjin 300072, China ∥ School of Chemistry and Chemical Engineering, University of the Chinese Academy of Sciences, Beijing 100049, China ‡

ABSTRACT: Most porous media (just like catalyst pellets) have complicated pore structures, and understanding the coupling of the diffusion and reaction processes in these pores is very important for improving their performance. In this work, a diffusion factor (D) and a reaction factor (R) are proposed to quantitatively describe the diffusion and reaction performance in these pores respectively at molecular level. The yield in unit time is used to quantify their productivity and is expressed as the product of D and R. Molecular dynamic simulations with the hard-sphere algorithm are carried out to study the reaction-diffusion coupling in several simple pore structures with the same volume, such as straight, T-shaped, and cross-shaped pores. The reaction formula based on activation energy is given for a simple irreversible reaction process from A to B. In terms of the proposed factors, D and R, analysis on the simulation results shows clearly that the overall productivity of these pore structures depends on the competition of D and R, which are both determined by the size and shape of the pore structures. The results demonstrate the effectiveness of the simulation approach used for evaluating the performance of the simple pore structures for simple reactions and the potential of its application in more complicated and practical cases. It also suggests the effectiveness of the proposed factors, D and R, for charactering the diffusion and reaction processes at molecular level.

1. INTRODUCTION Catalysts have wide applications and vital importance in many chemical processes, including fluid catalytic cracking (FCC),1−3 petroleum refining,4 polymerization,5 fuel cells6,7 and environmental protection.8 It is recognized that the performance of a single catalyst pellet (or particle) is significantly affected by the processes of diffusion of multicomponent, adsorption of reactant, reaction at surface, and desorption of reactant and product from the catalyst surface.9 In practice, the majority of catalyst is synthesized by the dispersion of nanoparticles containing active components on the nanoporous catalyst support. In general, a catalyst pellet with larger inner surface will contain more active sites and be beneficial to the catalytic performance. Obviously, for a given porosity, the smaller pores a catalyst pellet contains, the larger its inner surface area. However, too many small pores will cause diffusion difficulties. With the considerable progress in theoretic and experimental research, it is an increasing demand to develop new and efficient catalyst based on the understanding about the catalytic performance at the molecular level. Concerning the design of catalysts, one focus is on the design of catalytic surface and active sites, the other is on the design of catalyst pellets and their porous structures.10 Actually, the catalysts employed in © XXXX American Chemical Society

industry have very complicated multiscale porous structures, with pores of different sizes (micro-, meso- and macro-pores) and shapes coexisting. However, how to quantitatively describe the relationship between the reaction and diffusion process in restricted space, such as catalyst pellets, micro/nano-reactors and membrane reactors, is unclear yet, despite great efforts to optimize the pore networks.11 In practice, many factors need to be taken into account, such as porosity, specific surface area, connectivity, size distribution and geometrical shapes of the pore. In the past few years, complicated models of pore structures have been proposed, such as random oriented model,12 Bethe lattices,13,14 dual site-bond lattice model15 and so on. For example, a threedimensional crystallite-pore network model16 was established to study the multicomponent diffusion and reaction in microporous media with coexisting intracrystalline micro- and Special Issue: Tribute to Keith Gubbins, Pioneer in the Theory of Liquids Received: July 21, 2017 Revised: September 8, 2017

A

DOI: 10.1021/acs.langmuir.7b02559 Langmuir XXXX, XXX, XXX−XXX

Article

Langmuir

2. COMPUTATIONAL MODEL AND METHOD 2.1. Molecular Simulation Method. For HS model, event-driven algorithms28 have shown much higher efficiency for the diffusion of dilute gas in porous media than traditional time-driven algorithms if the number of molecules involved is not too large,29 and is hence employed in this study. In this approach, molecules are regarded as spheres with step potential that collide with each other instantly. When two hard spheres with mass (m1, m2), position (P1, P2), and velocity (V1, V2) collide after a free flight, their velocities after collision (V1′, V2′) are updated as

meso-pores. In their work, the transport was described by the Maxwell−Stefan surface diffusion model and the Knudsen diffusion model, respectively. And the reaction rate was calculated based on the coverage of the xylene isomers. Cao et al.17 calculated the diffusion coefficient by solving gas diffusion equations to study the gas diffusion in fractal porous. Depending on the unified Maxwell−Stefan diffusion theory, Li et al.18 proposed a multiregion model to study the reactiondiffusion processes in porous catalyst pellets, taking into account the size and spatial distribution of crystal particles in the pellet. Johannessen et al.19 introduced a continuum model of distributor channels network in nanoporous catalyst pellet to investigate its optimal structure. More detailed description of three-dimensional random network models and other models of catalyst support pore structure can be found in reviews by Keil.20,21 However, in these studies the coupling between diffusion and reaction and their effect on the overall productivity have not been the focus of interest. To study the reaction-diffusion processes within the complicated porous media, previous models have mostly based on dimensionless groups, such as Damköhler number, Thiele modulus, and the corresponding external/internal effectiveness factors,22−25 which were defined based on Fick’s diffusion law using empirical effective diffusion coefficients and idealized pellet shapes. As gas diffusion in complicated pores does not follow Fick’s law,17 and actual pore microstructures are not considered directly,23,24 these models and indexes, though helpful, are often used on a phenomenological and qualitative basis.26,27 Since the pore structure in actual porous media is too complicated to be described integrally and elaborately, it is meaningful to first focus on some typical structures at microscale and try to characterize the coupling of reaction and diffusion in a more quantitative way. For this purpose, molecular dynamics (MD) simulation can be used as a powerful tool in that the complicated diffusion and reaction processes are reproduced by the basic motion of the microscopic elements, rather than numerically computed from their material properties or statistical correlations, which are hard to obtain at this scale. However, as a large number of configurations and conditions have to be investigated for a quantitative understanding, reasonable models and efficient algorithms are both important for such a study. In this work, we target on a modest goal along this direction, that is, for a given volume of pore structure, we try to answer which geometry can produce the highest yield under given conditions, or in other words, which geometry is the most efficient. Though porous catalyst particles are taken as one of the background of this study, it is not limited to that. As it is impossible to exhaust all possible geometries in MD simulations, some typical simple shapes are investigated with a wide range of parameters first, which includes straight, Tshaped and cross-shaped pores. To be more concentrated on the primary factors for the coupling between reaction and diffusion, rather than the reaction itself, an extremely simplified reaction kinetics is designed and simple hard-sphere (HS) molecules are employed as reactant. This study is, on the one hand, in order to examine whether such simplifications still contain enough physics to obtain reasonable results, and on the other hand, to explore how diffusion and reaction coupling can be characterized by quantitative factors defined at microscales.

⎧ 2m2 (V1 − V2)·(P1 − P2) (P1 − P2) ⎪ V1′ = V1 − m1 + m2 |P1 − P2|2 ⎪ ⎨ ⎪ 2m1 (V1 − V2)·(P1 − P2) (P1 − P2) ⎪ V2′ = V2 + m1 + m2 |P1 − P2|2 ⎩

(1)

Many light gas molecules can be reasonably simplified as HS under normal and dilute conditions where an effective HS diameter σHS for the diffusion process can be calculated according to Ben-Amotz et al.,30 that is −1/6 ⎡ ⎛ T * ⎞1/2 ⎤ σHS(T *) ⎟ ⎥ = 1.1532⎢1 + ⎜ ⎝ 0.527 ⎠ ⎥⎦ ⎢⎣ σLJ

(2)

where, T* = kBT/εLJ, kB, and T represent the Boltzmann constant and the absolute temperature, respectively. The parameters σLJ and εLJ are from Lennard-Jones (L-J) potential, and they are the finite distance at which the interparticle potential is zero and the depth of the potential well, respectively. That means, it can still be used as a reasonable model for real molecules in gaseous state (and under simplified reaction scenarios). In the event-driven algorithm for HS, the elastic collisions are treated analytically, and all numerical errors, except the rounding error, can be avoided. Hence, the energies and momentums of the molecules could be preserved to machine accuracy. The HS approach has been widely used in the research of microfluidics. In our previous work,28,31 HS model has been applied to study the diffusion of gas in porous structures with reasonable results, which can be served as a validation of this method. More details of the HS approach can be found in these earlier publications. 2.2. Diffusion and Reaction Model. For a typical gas− solid chemical reaction, the reactant experiences a series of processes after entering the pore, including diffusion to the surface, adsorption of reactant, collision with the surface, reaction at surface, desorption from the surface and diffusion out of the pore. Sometimes, the reactant could pass through the pore without interacting with the surface. For simplicity, only two processes, diffusion and reaction, are considered in our model. Meanwhile, as we do not pay attention to a specific chemical reaction, a simple first-order irreversible isothermal reaction process from A to B (A → B) was considered in the present work. The active sites are assumed to be uniformly distributed over the entire surface of the pore. Based on the concept of activation energy, only collisions between molecules with high enough energy will lead to reactions when reactant molecules collide with the pore surface. Since there is no potential energy in the simulation model, the average energy of 1 reactant molecules can be expressed by 2 mv0 2 , where v0 is the thermal velocity obtained from the given temperature T. If the B


Article

Langmuir 1

(L), width (W), and height (H). The branches of the T-shaped and cross-shaped structures are described by their length (l), width (w) and height (h). If not specified, the main pore is open at its two ends and the branches are closed. Therefore, the reactant and product could diffuse inside the pore between the two ends of main pores. As shown in Figure 1, periodic boundary conditions are applied at the two ends, so once a molecule leaves one end, it immediately enters the other end. However, if it is a product (B) molecule, it will enter as a reactant (A) molecule. The pore volume (V = 9047153), molecule number (N = 360000) and energy (E) are the same for all the structures and remains unchanged. As potential energy is absent for HS molecules, the pressure (P) and temperature (T) were also conserved under the given isothermal reaction. At the initial state of simulation, all molecules in the pore are reactant and distributed uniformly with a fully random initial

energy of a reactant molecule ( 2 mvt 2 ) is higher than a certain 1

value, i.e., 2 mv0 2 , the reaction is considered to occur. That is, when a reactant (A) molecule collides with the pore surface, if its velocity vt satisfies vt 2 ≥

1 2 v0 η

(3)

(where η is the reaction coefficient), it will turn into a product (B) molecule. Obviously, larger η is easier for the reaction to take place. In fact, the reaction activation energy could be altered by changing η. 2.3. Characterizing Reaction and Diffusion Performance. The Thiele modulus proposed at macro-scale can be understood as the dimensionless ratio of the reaction and diffusion rates.32 The diffusion performance refers to the ability of the pore to provide reactant to react on the surface. The reaction performance refers to the ability of the surface to consume the reactant without mass transfer limitations. Following this understanding, the diffusion performance of a pore structure at microscale may be characterized by a diffusion factor (D) defined as the number of reactant colliding with unit pore surface area in unit time (just collisions, not necessarily reaction). In the same spirit, the reaction performance of the pore structure can be characterized by a reaction factor (R) with

R = P × Sin

velocity of v0 =

, which are then relaxed to equilibrium.

In addition, each molecule has a flag to indicate whether it is a reactant (A) or a product (B) molecule. Although the simulations can be conducted without referring to any specific system, and a dimensionless analysis can provide more general results, it is also helpful to understand its practical relevance by mapping the dimensionless parameters to those of some familiar systems. In this way, the simulated systems in this study can be approximated by the gas methane at 723 K and 9.98 MPa with σLJ = 1.0918 (3.7275 Å) and εLJ/kB = 148.99 K according to ref 33. The representative Knudsen number (for a cubic pore with the given volume) in this work is about 0.028, so the pore scale studied is mainly in the bulk and Knudsen diffusion regime. The simulation time step is 0.03 and the total simulation time is no less than 21753 to ensure that the system reaches a steady state. The relevant results are stored at every 100th step.

(4)

where P is the probability (or ability) of reaction at the active sites when collision occurs, which is obviously related to the reaction coefficient η and Sin is the total pore surface area. Finally, the overall performance of the pore structure, or its productivity, can be characterized by the number of product molecules diffusing out of the pore in unit time, which is defined as the yield (or product flux) Y. Obviously, Y=D×R

3kBT m

(5)

3. RESULTS AND DISCUSSIONS 3.1. Properties of the Steady State. In the simulations, the diffusion factor D, the reaction factor R and the yield Y are obtained in the steady state. As indicated in Figure 2, to obtain better statistical results, the number of reactant/product

Based on these definitions, the coupling between reaction and diffusion may be investigated quantitatively. To further facilitate the discussion and understanding, these two factors and other parameters hereafter are nondimensionalized by the molecular mass m and diameter σHS, and characteristic time σHS/v0 for mass, length and time, respectively. 2.4. Simulation Setup. The pore structures in this work are depicted in Figure 1. The main pore of the straight, Tshaped and cross-shaped structure is described by its length

Figure 2. Number of the reactant molecules reaching the surface (Nd), the reactant molecules reacting at the surface (Nr) and the product molecules diffusing out of the pore (Ny) in unit time versus the simulation time in straight pore (H = W = 176, L = 293, η = 0.25).

Figure 1. Schematic diagram of the pore structures. C


Article

Langmuir

in slice C, respectively, where the densities are averaged over the previous 20000 steps for each snapshot, the reactant near the surface gradually decreases (Figure 4) but the product gradually increases (Figure 5) as the simulation goes. Under the steady state, their distributions also remain unchanged statistically, which is the result of the coupling between the reaction and diffusion processes. The densities of the reactant and product in slices A, B, and D at t = 21753 are then shown in Figures 6 and 7. It can be found that most reactant appears near the two ends of the pore (see slice B in Figure 6). In the cross-section at the end (see slice D in Figure 6) or in the middle (see slice A in Figure 6), the reactant tends to appear at the central region of the section because they are easier to react near the surface. On the contrary, most product appears in the middle region of the pore (see slice B in Figure 7). In the cross-section (see slice A and D in Figure 7), the product tends to appear near the surface, especially at the corner. These results are certainly consistent with our expectation under the given reaction rules. 3.2.2. Effect of Size on the Coupling of Reaction and Diffusion. As the pore volume V remains the same for different pore structures, according to the formula Y = D × R, the change of the pore geometry will lead to different performance of the reaction and diffusion processes. It is hence interesting to know what kind of geometry can achieve better performance. 3.2.2.1. Effect of Pore Length with a Square Cross-Section (H = W). In this section, the reaction coefficient η is kept at 0.25 and the cross-sectional area Sc = H × W is varied together with the length of the pore (L) under the constraint of H = W. As

molecules are sampled at every 10000 steps and then averaged. It is found that Nd, Nr, and Ny fluctuated around their constant values after about t = 12430, which indicates that the simulated system reaches the steady state. The time needed to reach the steady state is related to the reaction coefficient η, the pore size, and structure. 3.2. Straight Pore. 3.2.1. Density Distribution. Variation of the number density distribution of the reactant/product is a direct indication of the reaction and diffusion processes, for which several slices inside the straight pore are selected for analysis, as shown in Figure 3. The thickness of each slice is

Figure 3. Diagram of the sampled slices in the straight pore (slices A, B, and C are in the middle position along the a-, b- and c-directions, respectively, and slice D is located at the open end along the adirection).

about 11.72, with the sampling mesh sized at 5.86 × 5.86 × 11.72. As shown in Figures 4 and 5 for the reactant and product

Figure 4. Number density distribution of the reactant in slice C at different simulation time (t) for straight pore (H = W = 176, L = 293, η = 0.25). D


Article

Langmuir

Figure 5. Number density distribution of the product in slice C at different simulation time (t) for straight pore (H = W = 176, L = 293, η = 0.25).

Figure 6. Number density distribution of the reactant at simulation time t = 21753 (steady state) in different slices of the straight pore (H = W = 176, L = 293, η = 0.25).

Figure 7. Number density distribution of the product at simulation time t = 21753 (steady state) in different slices of the straight pore (H = W = 176, L = 293, η = 0.25).

E


Article

Langmuir shown in Figure 8, with the increase of L, Y first increases and then decreases, reaching its maximum at L = 293 roughly. To

increases, it is difficult to diffuse out. As shown in Figure 10c, the product almost fills the whole pore. Therefore, the optimal coupling can be achieved only when the reaction and diffusion reach a compromise in the pore at a specific length (Figure 10b). 3.2.2.2. Effect of the Cross-Sectional Shape. The geometrical structure of the pore can affect the reaction and diffusion processes. An interesting question is how the crosssectional shape affects the overall performance when the pore length remains constant. For this purpose, η is set to 0.25 as before and the optimal pore length for the square cross-section, L = 293, is kept and hence the cross sectional area Sc = 30887 remains constant. The ratio of the pore length to width of the cross-section λH/W = H/W is varied from 1 to 3600, which means the narrow aspect is only 3 molecular diameters in size and its further decrease is impractical. As the cross-section changes from square to narrow rectangle, Y increases monotonically as shown in Figure 11. That is to say, the pore with a square cross-section has the lowest performance. To understand this result, the variation of the diffusion factor D and the reaction factor R versus λH/W is shown in Figure 12. It is easy to know that the surface area increases with the increase of λH/W and R. However, as the width of the pore reduces, the molecular motion in this direction is limited, so D decreases with the increase of λH/W. These opposite variation tendencies indicate the compromise between the diffusion and the reaction processes. Since Y is increasing, we can infer that the influence of λH/W on R exceeds that on D. It should be pointed out that when the pore becomes too narrow, the diffusion mechanism may change and the adsorption and desorption may dominate, which are not considered in the present study. 3.2.3. Effects of Changing the Reaction Coefficient η. As mentioned before, the effects of activation energy can be investigated by adjusting the reaction coefficient η. To study this effect, the straight pore with the same given volume (H = W = 176, L = 293) is used and its length still corresponds to the maximum value of yield Y in Figure 8. As shown in Figure 13, Y increased rapidly in the range of η < 1.0, and then slowly when η > 1.0. As larger η indicates lower activation energy and more reactions occurring in unit time, R increases with the increase of η, which is verified in Figure 14. However, Y cannot always increases rapidly because both the velocity and density distributions of the reactant are bounded. Since the pore structure has not changed, the physical diffusion process will not be affected. However, the diffusion factor D decreased rapidly when η < 0.5, and then decreases slowly with the further increase of η, as shown in Figure 14. The reason is that when the activation energy gets lower, more reactant will be converted into product, which results in a marked reduction in the number of reactant. So the number of collisions between the reactant and the surface decreases. Because D was defined as the number of reactant colliding with the unit surface in unit time, D decreases with the increase of η. That means that the diffusion coefficient as a material property is not a comprehensive index of the diffusion effect, as it is not coupled to the reaction process, while the diffusion factor D defined in this study can take a more effective role. 3.2.4. Effects of Opening of the Pore. In the catalyst pellet, the straight pores are not always open at both ends. In this section, we will discuss the coupling of reaction and diffusion in three different types of pores, which all have square cross sections as shown in Figure 15. The first one (N-type) is open

Figure 8. Variation of the yield Y versus the length L of the straight pore (H = W, η = 0.25).

better understand this point, the variations of the diffusion factor (D) and the reaction factor (R) versus L are shown in Figure 9. As L increases, D decreases, and R increases gradually.

Figure 9. Variation of the diffusion factor D and the reaction factor R versus the length L of straight pore (H = W, η = 0.25).

For a given volume, the increase of L leads to the increase of pore surface area and hence the active sites, therefore, the increase of R is reasonable in view of its definition. On the other hand, the increase of L means the decrease of the crosssectional area, which is not beneficial to the diffusion of reactant and product. Although more product can be produced, they are difficult to diffuse out of the pore, and then will occupy the active sites and even hinder the diffusion of the reactant to the active sites. Therefore, the diffusion factor D apparently reduces with the increase of L. In addition, the number density distributions of the product in slice C with different L under the steady state are shown in Figure 10. When L is small, the product is easy to diffuse out of the pore because of the relatively large cross-sectional area, but less product is produced because the surface area is smaller. Therefore, the product is relatively less and distributed near the surface (Figure 10a). However, when L is larger, although the product generated F


Article

Langmuir

Figure 10. Number density distribution of the product at simulation time t = 21753 in slice C with different L under the steady state (H = W, η = 0.25).

Figure 13. Variation of yield Y versus the reaction coefficient η (H = W = 176, L = 293).

Figure 11. Variation of Y versus λH/W (L = 293, η = 0.25).

at both ends, which is the one used above, the second one (Etype) is closed at one end, and the third one (M-type) has a wall in its middle. It is noted that the newly added wall is regarded as the active inner surface, on which reactions can also take place. For these three types of pores, the variation of the yield Y versus the length L is investigated and the results are summarized in Figure 16. As L increases, the Y values of the E-type and the M-type pores decrease monotonically, but the Y of the N-type increases first and then decreases. When L is smaller (about L < 164), the performance of the M-type pore is the best, but the performance of the N-type pore is the worst. Figure 16 also shows that when L > 659 roughly, the performance of the M-type and the N-type pores is almost the same and better than that of the E-type. When L is very large (about L > 10000), the three types of pores will have almost the same performance.

Figure 12. Variation of D and R versus λH/W (L = 293, η = 0.25).

G


Article

Langmuir

Figure 14. Variation of diffusion factor D and reaction factor R versus the reaction coefficient η (H = W = 176, L = 293). Figure 16. Variation of the yield Y versus the length L of the three types of straight pores (η = 0.25).

To better understand the performance of the three types of pores, the diffusion factor D and the reaction factor R are shown in the Figure 17. It is shown that as L increases, D of the three types of pores decrease. The R deceases first and then increases for the M-type and the E-type, but increases monotonically for the N-type. When L > 659, the R of the three types of pores have almost the same value. Meanwhile, D of the M-type and the N-type is almost the same and larger than that of the E-type, which is in good agreement with the variation of Y in Figure 16. That means, long and narrow (about L > 659) pores are diffusion dominated while shorter and wider pores are reaction dominated (about L < 164). When the length is within the range of 164 and 659, both diffusion and reaction are dominant factors, which results in complex coupled performance. In addition, the fact that the M-type pore has the best performance for shorter length and the N-type and M-type pores have the same diffusion performance for longer length indicate that most gas molecules do not enter the pore from one end and diffuse out from the other end. That is to say, for a catalyst pellet with a given porosity, the effective distance the gas molecules can penetrate into the pore is limited, which is consistent with the conclusion obtained by considering the reaction and diffusion using Thiele modulus.23 3.3. T-Shaped Pore. In addition to the straight pore, the Tshaped pore is also a commonly used simplified model, as shown in Figure 1. Since there are several structural parameters that can be optimized, here, only the effect of the height of branch (h) on the reaction-diffusion coupling performance is investigated. In the simulations, the volume of the T-shaped pore was the same as that of the straight pore. The cross-sectional area of the main pore in the T-shaped structure also corresponds to the maximum value of yield Y of the straight pore in Figure 8.

Figure 17. Variation of diffusion factor D and reaction factor R with the length L of three types of straight pores (η = 0.25).

The ends of the branches closed. When h changes, the length of the main pore (L) will change accordingly to ensure that the volume remains constant. Besides, the reaction coefficient (η = 0.25), the width (W = 176) and height (H = 176) of the main pore, and the width (w = 176) and length (l = 88) of the branches remain unchanged. As shown in Figure 18, with the increase of h, Y increases first and then decreases, reaching its maximum at h = 59 roughly. This is also the result of the compromise of the diffusion and reaction. As shown in Figure 19, R increased with the increase of h, which is mainly due to the increase of the surface area. In addition, the increase in h leads to a decrease in L. According to the simulation results of the straight pore, D should increase. However, D decreases with the increase of h in Figure 19. It can

Figure 15. Schematic diagram of the different structure types of straight pores with square sections. H


Article

Langmuir

however, the overall performance mainly depends on the size of the main pore. 3.4. Cross-Shaped Pore. The cross-shaped pore is also a simplified model which can be regarded as a basic unit of complex pore structure, as shown in Figure 1. First, similar to the study of the T-shaped pore above, the effect of the branch height (h) on the reaction-diffusion coupling is investigated. 3.4.1. Effects of the Branch Height. In the simulations, the volume of the cross-shaped pore was the same as that of the straight pore. The ends of the branch are closed. When h changes, the length of the main pore (L) will change accordingly to ensure that the volume remains constant. The dimensions of the two branches are always the same and the cross sectional area of the main pore in the cross-shaped structure is also the same as that of the straight pore with the maximum value of yield Y in Figure 8. Besides, other structural sizes remain unchanged, including the width and height (H = W = 176) of the main pore, the width (w = 176) and length (l = 29) of the two branches, as well as the reaction coefficient (η = 0.25). The effect of branch height h on the yield Y is similar to that in the T-shaped pore and Y reached its maximum value when h is about 88 (see Figure 21). The compromise of the diffusion

Figure 18. Variation of the yield Y versus the height (h) of branch in the T-shaped pore (η = 0.25, H = W = w = 176, l = 88).

Figure 19. Variation of the diffusion factor D and the reaction factor R versus the height (h) of branch in T-shaped pore (η = 0.25, H = W = w = 176, l = 88).

be explained that the increase of the branch height blocks the diffusion process. As shown in Figure 20, the product in the branch is difficult to diffuse out of the main pore, which results in lower overall performance. Besides, the relative variation of Y shown in Figure 18 is small. This is because the variation of the branch height leads to slight variation of the main pore length,

Figure 21. Variation of the yield Y versus the branch height (h) in the cross-shaped pore (η = 0.25, H = W = w = 176, l = 29).

and reaction when changing h is shown in Figure 22. Similar to the case of T-shaped pore, the product in the two branches is

Figure 20. Number density distribution of the product at simulation time t = 21753 in Slice C in the T-shaped pore with different branch heights (h = 12, 47 and 117, respectively) under the steady state (η = 0.25, H = W = w = 176, l = 88). I


Article

Langmuir

Figure 22. Variation of the diffusion factor D and the reaction factor R versus the branch height h in the cross-shaped pore (η = 0.25, H = W = w = 176, l = 29).

Figure 23. Number density distribution of the product at simulation time t = 21 753 in slice C of the cross-shaped pore with different branch heights (h = 18, 88 and 193, respectively) under steady state (η = 0.25, H = W = w = 176, l = 29).

difficult to diffuse out of the main pore, as shown in Figure 23. However, it is also found that the maximum value of Y obtained in a cross-shaped pore is larger than that in a T-shaped pore, and that in a straight pore is the smallest. So, the appearance of branches can reduce the diffusion performance, but increase the reaction performance due to the increase of the reactive surface, and the overall performance becomes better. 3.4.2. Effects of Branch Numbers. From the simulation results above, it can be found that adding branches on a similar main pore structure can improve the performance. To explore whether even more branches can improve the performance further, two branches are added at a time, and all branches are of the same structural size. For convenience, these structures are denoted as cross-2, cross-4, and cross-6, respectively. Similar to the previous optimization process, the optimal performance is obtained by changing the branch height (h) and the length (L) of the main pore. Each maximum value of Y obtained from different structure corresponds to a different h. Then the maximum values of Y from different structures are shown in Figure 24. It is shown that with more branches, the overall performance is higher. The corresponding diffusion factor D and the reaction factor R are shown in Figure 25. More

Figure 24. Best yield Y in straight, T-shaped, and different crossshaped pores (η = 0.25, H = W = w = 176).

branches will reduce the diffusion performance but increase the reaction performance. Therefore, for a catalyst pellet with a J


Article

Langmuir

Figure 25. Variation of the diffusion factor D and the reaction factor R in straight, T-shaped, and cross-shaped pores (η = 0.25, H = W = w = 176).

Figure 26. Variation of the yield Y versus the branch height h in the cross-shaped pore with open branches (η = 0.25, H = W = w = 176).

given porosity, more branches can improve its performance, but its structure must be optimized. Some relevant studies using continuum-based reactiondiffusion equations can also be found for the simple pore structures similar to our study above. For example, Johannessen et al.,19 Gheorghiu and Coppens34 and Wang et al.35 have studied the optimal networks of distributor channels in porous catalysts. In their work, an optimal structure (maximum yield) can also be obtained in a catalyst pellet of a given geometry and volume. Their work is significant for the design of catalyst structure and will help in improving the performance of catalyst. However, it is not possible to make a direct comparison between the results obtained by HS simulation and by the continuum model because only a simple virtual reaction rather than a specific reaction is considered in the present work, and it is difficult to give suitable parameters for continuum model. In our future work, the simulation model and computing capability will be further improved, then direct comparison will become possible. 3.4.3. Effects of Open Branches in Cross-Shaped Pores. The pores in the catalyst pellets are actually interconnected, forming a pore network structure. The cross-shaped pore structure with open branches is more realistic for studying actual performance. In the simulations, the cross-shaped pore and the structural optimization process are same to the above. The only difference is that the ends of the branches are now open. The variation of the yield Y versus the branch height h is shown in Figure 26. It is found that the maximum Y in the structures with open branches is about 1.3 times of those with closed ends, as shown in Figure 27. Meanwhile, the optimal height (h) of the open branches is 2.6 times of the closed branches.

Figure 27. Yield Y in cross-shaped pores with different number of branches (η = 0.25, H = W = w = 176).

by the presence or absence of an inner wall in its middle, due to the limited effective diffusion distance of gas molecules. An optimal pore length exists for the same pore volume under given conditions. Besides, larger aspect ratio of the pore cross sections seems to improve the overall performance of the straight pores. The results also suggest that the performance of cross-shaped pores is better than that of T-shaped and straight pore with the same pore volume. Moreover, a main pore with more open branches will have even better performance, but their structures must be optimized to achieve maximum performance. The results demonstrate the effectiveness of the simulation approach used for evaluating the performance of the simple pore structures for simple reactions and the potential of its application in more complicated and practical cases. It also suggests the effectiveness of the proposed factors, D and R, for characterizing the diffusion and reaction processes at a molecular level. However, the present work is based on very simple pore structures and reaction kinetics, and surface diffusion and adsorption/desorption are not considered yet. The effects of

4. CONCLUSIONS Molecular dynamic simulations with hard-sphere model are carried out to investigate the reaction-diffusion coupling in several simple pore structures with the same volume, including straight, T-shaped, and cross-shaped pores. The simulation results demonstrated that the overall performance of the pore structures depends on the competition of the diffusion and reaction factors. It is found that the diffusion factor of a sufficiently long and narrow straight pore is almost unaffected K


Article

Langmuir λH/W

Knudsen number and/or molecular concentration are also not investigated. A large amount of future work is, therefore, needed, with which we hope some dimensionless criterion can be proposed based on the D and R defined in this work to better understand the coupling between reaction and diffusion processes.

■

■

ratio of the pore length to width of the cross-section

REFERENCES

(1) Errekatxo, A.; Ibarra, A.; Gutierrez, A.; Bilbao, J.; Arandes, J. M.; Castaño, P. Catalytic Deactivation Pathways During the Cracking of Glycerol and Glycerol/VGO Blends under FCC Unit Conditions. Chem. Eng. J. 2017, 307, 955−965. (2) Mance, D.; van der Zwan, J.; Velthoen, M. E.; Meirer, F.; Weckhuysen, B. M.; Baldus, M.; Vogt, E. T. A DNP-Supported SolidState NMR Study of Carbon Species in Fluid Catalytic Cracking Catalysts. Chem. Commun. 2017, 53 (28), 3933−3936. (3) Vogt, E. T. C.; Weckhuysen, B. M. Fluid Catalytic Cracking: Recent Developments on the Grand Old Lady of Zeolite Catalysis. Chem. Soc. Rev. 2015, 44 (20), 7342−7370. (4) Abdalla, A.; Arudra, P.; Al-Khattaf, S. S. Catalytic Cracking of 1Butene to Propylene Using Modified H-ZSM-5 Catalyst: A Comparative Study of Surface Modification and Core-Shell Synthesis. Appl. Catal., A 2017, 533, 109−120. (5) Chen, M.; Chen, C. Rational Design of High-Performance Phosphine Sulfonate Nickel Catalysts for Ethylene Polymerization and Copolymerization with Polar Monomers. ACS Catal. 2017, 7 (2), 1308−1312. (6) Yu, W.; Tao, J.; Yu, X.; Zhao, S.; Tu, S.-T.; Liu, H. A Microreactor with Superhydrophobic Pt−Al2O3 Catalyst Coating Concerning Oxidation of Hydrogen Off-Gas from Fuel Cell. Appl. Energy 2017, 185, 1233−1244. (7) Antolini, E. Alloy vs. Intermetallic Compounds: Effect of the Ordering on the Electrocatalytic Activity for Oxygen Reduction and the Stability of Low Temperature Fuel Cell Catalysts. Appl. Catal., B 2017, 217, 201−213. (8) Ibarra, Á .; Veloso, A.; Bilbao, J.; Arandes, J. M.; Castaño, P. Dual Coke Deactivation Pathways During the Catalytic Cracking of Raw Bio-Oil and Vacuum Gasoil in FCC Conditions. Appl. Catal., B 2016, 182, 336−346. (9) Keil, F. J. Multiscale Modelling in Computational Heterogeneous Catalysis. Top. Curr. Chem. 2011, 307, 69−107. (10) Keil, F. J.; Rieckmann, C. Optimization of Three-Dimensional Catalyst Pore Structures. Chem. Eng. Sci. 1994, 49, 4811−4822. (11) Ye, G. H.; Zhou, X. G.; Zhou, J. H.; Yuan, W. K.; Coppens, M.O. Influence of Catalyst Pore Network Structure on the Hysteresis of Multiphase Reactions. AIChE J. 2017, 63 (1), 78−86. (12) Johnson, M. F. L.; Stewart, W. E. Pore Structure and Gaseous Diffusion in Solid Catalysts. J. Catal. 1965, 4 (2), 248−252. (13) Beeckman, J. W.; Froment, G. F. Catalyst Deactivation by Site Coverage and Pore Blockage: Finite Rate of Growth of the Carbonaceous Deposit. Chem. Eng. Sci. 1980, 35 (4), 805−815. (14) Reyes, S.; Jensen, K. F. Estimation of Effective Transport Coefficients in Porous Solids Based on Percolation Concepts. Chem. Eng. Sci. 1985, 40 (9), 1723−1734. (15) Armatas, G. S. Determination of the Effects of the Pore Size Distribution and Pore Connectivity Distribution on the Pore Tortuosity and Diffusive Transport in Model Porous Networks. Chem. Eng. Sci. 2006, 61 (14), 4662−4675. (16) Ding, W.; Li, H.; Pfeifer, P.; Dittmeyer, R. Crystallite-Pore Network Model of Transport and Reaction of Multicomponent Gas Mixtures in Polycrystalline Microporous Media. Chem. Eng. J. 2014, 254, 545−558. (17) Cao, L.; He, R. Gas Diffusion in Fractal Porous Media. Combust. Sci. Technol. 2010, 182 (7), 822−841. (18) Li, H.; Ye, M.; Liu, Z. M. A Multi-Region Model for Reaction− Diffusion Process within a Porous Catalyst Pellet. Chem. Eng. Sci. 2016, 147, 1−12. (19) Johannessen, E.; Wang, G.; Coppens, M.-O. Optimal Distributor Networks in Porous Catalyst Pellets. I. Molecular Diffusion. Ind. Eng. Chem. Res. 2007, 46 (12), 4245−4256. (20) Keil, F. J. Diffusion and Reaction in Porous Networks. Catal. Today 1999, 53 (2), 245−258. (21) Keil, F. J. Complexities in Modeling of Heterogeneous Catalytic Reactions. Comput. Math. Appl. 2013, 65 (10), 1674−1697.

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected] (C.X. Li). *E-mail: [email protected] (W. Ge). ORCID

Wei Ge: 0000-0002-4997-5736 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors gratefully acknowledge financial support from the National Natural Science Foundation of China grant Nos. 2146232, 21225628, 91434201, 91434102, and Key Research Program of Frontier Sciences, Chinese Academy of Sciences grant No. QYZDJ-SSW-JSC029, and the grants from the Center for Mesoscience, Institute of Process Engineering, Chinese Academy of Sciences, COM2015A006 and from CAS Interdisciplinary Innovation Team.

■

LIST OF SYMBOLS D diffusion factor R reaction factor Y yield m1, m2 mass σHS effective diameter of hard sphere kB Boltzmann constant T absolute temperature σLJ finite distance at which the interparticle potential is zero in Lennard-Jones (L-J) potential εLJ depth of the potential well in Lennard-Jones (L-J) potential vt velocity at collision v0 thermal velocity η reaction coefficient P probability (or ability) of reaction at each active site when collision occurs Sin total inner surface area of the pore L length of main pore W width of main pore H height of main pore l length of branches w width of branches h height of branches V pore volume N molecule number E energy P pressure Nd number of the reactant molecules reaching the surface in unit time Nr number of the reactant molecules reacting at the surface in unit time Ny number of the product molecules diffusing out of the pore in unit time t nondimensionalized simulation time Sc cross-sectional area L


Article

Langmuir (22) Karakaya, C.; Weddle, P. J.; Blasi, J. M.; Diercks, D. R.; Kee, R. J. Modeling Reaction−Diffusion Processes within Catalyst Washcoats: I. Microscale Processes Based on Three-Dimensional Reconstructions. Chem. Eng. Sci. 2016, 145, 299−307. (23) Blasi, J. M.; Weddle, P. J.; Karakaya, C.; Diercks, D. R.; Kee, R. J. Modeling Reaction−Diffusion Processes within Catalyst Washcoats: II. Macroscale Processes Informed by Microscale Simulations. Chem. Eng. Sci. 2016, 145, 308−316. (24) Förtsch, D.; Essenhigh, R. H.; Schnell, U.; Hein, K. R. G. On the Application of the Thiele/Zeldovich Analysis to Porous Carbon Combustion. Energy Fuels 2003, 17 (4), 901−906. (25) Thiele, E. W. Relation Between Catalyst Activity and Size of Particle. Ind. Eng. Chem. 1939, 31 (1), 916−920. (26) Essenhigh, R. H.; Klimesh, H. E.; Förtsch, D. Combustion Characteristics of Carbon: Dependence of the Zone I−Zone II Transition Temperature (Tc) on Particle Radius. Energy Fuels 1999, 13 (4), 826−831. (27) Essenhigh, R. H.; Förtsch, D.; Klimesh, H. E. Combustion Characteristics of Carbon: Influence of the Zone I−Zone II Transition on Burn-Out in Pulverized Coal Flames. Energy Fuels 1999, 13 (5), 955−960. (28) Zhang, C. L.; Shen, G. F.; Li, C. X.; Ge, W.; Li, J. H. HardSphere/Pseudo-Particle Modelling (HS-PPM) for Efficient and Scalable Molecular Simulation of Dilute Gaseous Flow and Transport. Mol. Simul. 2016, 42 (14), 1171−1182. (29) Akkaya, V. R.; Kandemir, I. Event-Driven Molecular Dynamics Simulation of Hard-Sphere Gas Flows in Microchannels. Math. Probl. Eng. 2015, 2015 (6), 1−12. (30) Benamotz, D.; Herschbach, D. R. Estimation of Effective Diameters for Molecular Fluids. J. Phys. Chem. 1990, 94 (3), 1038− 1047. (31) Li, Y. P.; Zhang, C. L.; Li, C. X.; Liu, Z. C.; Ge, W. Simulation of the Effect of Coke Deposition on the Diffusion of Methane in Zeolite ZSM-5. Chem. Eng. J. 2017, 320, 458−467. (32) Aris, R. Communication. Normalization for the Thiele Modulus. Ind. Eng. Chem. Fundam. 1965, 4 (2), 227−229. (33) Vrabec, J.; Fischer, J. Vapor-Liquid Equilibria of Binary Mixtures Containing Methane, Ethane, and Carbon Dioxide from Molecular Simulation. Int. J. Thermophys. 1996, 17 (4), 889−908. (34) Gheorghiu, S.; Coppens, M.-O. Optimal Bimodal Pore Networks for Heterogeneous Catalysis. AIChE J. 2004, 50 (4), 812−820. (35) Wang, G.; Johannessen, E.; Kleijn, C. R.; de Leeuw, S. W.; Coppens, M.-O. Optimizing Transport in Nanostructured Catalysts: A Computational Study. Chem. Eng. Sci. 2007, 62 (18−20), 5110−5116.

M


Simulation Study on the Reaction-Diffusion ... - ACS Publications

Recommend Documents