Research Article  Open Access
Improving Accuracy of Coarse Grid Numerical Solution of SolidSolid Reactions by Taylor Series Expansion of the Reaction Term
Abstract
Exothermic solidsolid reactions lead to sharp reaction fronts that cannot be captured by coarse spatial mesh size numerical simulations that are often required for largescale simulations. We present a coarsescale formulation with high accuracy by using a Taylor series expansion of the reaction term. Results show that such expansion could adequately maintain the accuracy of finescale behavior of a constant pattern reaction front while using a smaller number of numerical grid cells. Results for a onedimensional solidsolid reacting system reveal reasonable computational time saving. The presented formulation improves our capabilities for conducting fast and accurate numerical simulations of industrialscale solidsolid reactions.
1. Introduction
Modeling of reactive flow has diverse applications in engineering and science. Applications include heavy oil recovery processes, combustion in porous media, ground water flow, and transport and reaction processes in biofilms. An important class of reactions is that of solidsolid reactions. Modeling of solidsolid reactions is of significant interest in various industrial operations, including oxidation of metallic and nonmetallic mixed powders, cement industry [1], ferrites manufacturing, solidstate polymerization [2, 3], ceramic manufacturing [4], catalyst preparation [5], and drug storage (for more details, see a comprehensive review by Tamhankar and Doraiswamy [6]). Numerous investigations on modeling reactive flow have been reported over the years that greatly improved the understanding of such systems, where the emphasis has primarily been on approximate analytical solutions for special cases, fine grid direct numerical simulation, and upscaling of reactiontransport from pore scale to continuum scale [7–18].
Accurate numerical simulation of solidsolid reactions is a challenging task due to the multiscale nature of the physical phenomena. Physical processes involved in solidsolid reactive systems include diffusive (heat and mass) and reactive processes. Reactions in porous media intrinsically often take place at the small scale, causing development of subdiffusivescale concentration and temperature gradients, while heat and mass diffusive processes have scales orders of magnitude larger than the reactions. Largescale simulation of such coupled processes is computationally expensive due to limitations in computational resources. Therefore, a formulation that captures the subdiffusive scale improves our capabilities for conducting industrialscale simulation of the involved processes with high accuracy. In this paper, we provide a coarsescale formulation to capture the subgridscale phenomena appropriate for largescale numerical simulations. The paper is organized as follows. First, the finescale mathematical model used in this study is presented. Next, the coarsescale model is described. Then, application of the coarsescale model is given for a constant pattern reaction front, followed by summary and conclusions.
2. FineScale Conservation Equations
The exothermic solidsolid reaction is assumed to take place in a onedimensional semiinfinite domain. The system, which consists of a mixture of two solid materials, is initially at temperature . The concentration of one of the solids is in excess such that the reaction can be considered of first order with respect to the second solid. The initial concentration of the second solid is . At time , the initial temperature at is suddenly raised to ignition temperature high enough to initiate the reaction by local heating. The temperature dependence of the reaction rate is assumed to follow Arrhenius type behavior, and all physical properties are assumed to be constant. Mass diffusion is considered to be negligible, and the reaction front is planar and nonoscillatory. It is further assumed that the solid reacting mixture acts as an isotropic homogeneous system and that radiation effects are negligible [11]. The dimensionless energy and mass balances can be presented by parabolic coupled partial differential equations given by where the following dimensionless groups are used [11, 13]: where T is temperature, C is concentration of reactant, is density, is heat capacity, is the average thermal conductivity, k is the preexponential factor, E is activation energy, is the heat of reaction, R is the universal gas constant, t is time, is the scale temperature, and x is the spatial coordinate. We used a standard scaling available in the combustion literature to render the equations dimensionless [13, 19–21]. The scaling variable corresponds to an approximate measure of the heating zone length, and is a measure of the reaction front velocity [13].
The initial and boundary conditions are then given by The behavior of such a reacting convection diffusion system is primarily determined by the inverse of dimensionless activation energy and inverse of dimensionless heat of reaction, namely, and , respectively. Puszynski et al. [11] presented a detailed analysis of the frontal behavior of such a system. Small or values correspond to a reacting system with high activation energy or heat of reaction and vice versa. Large values of (small values of heat of reaction and/or small values of activation energy) lead to a degenerated combustion regime, whereas low values of (large values of heat of reaction and/or large values of activation energy) result in an oscillatory reaction front. Low values (large activation energy) and intermediate values result in a constant pattern profile regime. Choosing large grid cells results in significant smearing of the reaction front, as we will see in the following sections. In the subsequent part of the paper, we present the model formulation that is able to maintain the accuracy of fine grid behavior while using coarse grid cell sizes.
3. CoarseScale Conservation Equations
Numerical simulation of reactive front propagation with a large number of grid cells is computationally expensive and therefore we are interested in using a coarse model that preserves the accuracy of fine grid behavior. The coarsescale variable can be defined by where v represents the coarse grid.
We intend to represent the differential equations (2.1) with their corresponding coarse grid equivalents (3.2). In these equations, represents the “equivalent” reaction rate that would allow close agreement between the (2.1) set and the (3.2) set, correspondingly. The coarsescale numerical model can be expressed by [18] where are the finescale reaction rate and coarsescale average reaction rate, respectively. The reaction rate can be approximated by a Taylor series expansion around the coarsescale temperature and concentration. The above expansion can be represented in terms of deviation from the coarsescale variables as given by where and are the concentration and temperature deviations from the fine scale solution, respectively. The terms involving derivatives account for the information that is normally lost in a coarse grid numerical solution; those will be evaluated and accounted for herein. By using (3.1), the average coarsescale reaction rate can be expressed by where by definition the terms with first derivatives are dropped in the averaging process [18]. Equation (3.5) is a coarsescale representation of the reaction rate that includes both coarsescale concentration and temperature and their deviations from fine scale, namely, and , respectively. We do not have an explicit definition of the deviation terms and as functions of the coarsescale concentration and temperature. Instead, we intend to express the quantities and in a form proportional to and such that the final coarsescale model can be presented as a function of coarsescale temperature and concentration only. Similar approach has been used by Meile and Tuncay [18] for linear convectiondiffusionreaction problem. The common practice in a volume averaging method [14–17, 22] is to solve the closure problem, which is obtained by subtraction of the finescale and the coarsescale equations. Here, due to the complexity arising from the nonlinearity of the problem, subtraction does not lead to (an) equation(s) in terms of the perturbed quantities. However, consistent with the volume averaging method [14–16], an appropriate approximation for the deviation terms is linear proportionality of temperature and concentration deviations with their corresponding gradient. Using this approximation, one may write By substituting (3.9) in (3.5), the coarsescale equations (3.2) can be written as where is the proportionality constant or coarsescale parameter and is a function of coarse grid size. The coarsescale formulation is not closed so far, and the proportionality constant as a function of coarse grid size needs to be determined. Numerical experiments reveal that the term containing cross derivatives in the internal bracket in the righthand side in (3.10) is small as compared to the second term. Given that the cross derivatives are small as compared to the second term, it may be ignored. Therefore, the coarsescale formulation can be rendered as The coarse grid model represented by (3.11) is of the same form as (2.1) with the same dimensionless groups, and ; however, the reaction term is replaced by .
In order to find the proportionality constant as a function of coarse spatial mesh size, a series of numerical experiments needs to be conducted for different spatial mesh sizes for each specific reacting system (i.e., using fixed and ). The functionality ( versus spatial mesh size) for a specific reacting system can be determined by matching the coarse spatial mesh size solution with the fine scale or reference solution. In matching process, to estimate the deviations of the coarsescale model predictions from the finescale reference solution, we define a numerical error using the following expression as a measure of accuracy: where can be either temperature or concentration.
The proportionality constant is obtained by the minimization of the numerical error given by (3.12). By repeating the matching process for different number of spatial mesh sizes, the proportionality constant for each spatial mesh size can be obtained. Once the functionality of the proportionality constant with respect to coarse spatial mesh size is obtained, the coarsescale formulation (3.11) is complete and can be used for subsequent largescale simulation of a solidsolid reactive system. In the following, we use the above procedure for determining the parameter as a function of coarse spatial mesh size for some solidsolid reacting systems.
4. Application of the CoarseScale Formulation
Based on the scaling groups used in nondimensionalizing the governing equations, the frontal behavior of a solidsolid reaction depends on two parameters, namely, (inverse of dimensionless activation energy) and (inverse of dimensionless heat of reaction). The coarsescale formulation described previously is applied for two solidsolid reactions. The governing differential equations are discretized using an explicitintime finite difference approximation. A blockcentered scheme is used, where the diffusive flux is calculated based on grid block center values. Numerical simulations are conducted to determine the parameter appropriate for largescale numerical simulation of such reacting systems. The dimensionless constants for these reacting systems are given in Table 1. The data are taken from examples of solidsolid reactions given by Puszynski et al. [11]. For each reacting system, temperature at one end is rapidly increased to an ignition temperature. Since the temperature is scaled with the adiabatic temperature of the reaction , the minimum dimensionless temperature in the system is equal to the negative of the dimensionless heat of reaction . Figure 1 shows as a function of numerical spatial mesh size for different solidsolid reaction systems given in Table 1. In all cases, the conditions are selected such that a constant pattern profile exists. The condition for existence of a constant pattern reaction front can be predicted by the following expression [11]: where For , a reacting system demonstrates constant pattern reaction front propagation. The parameter for reactions studied is given in Table 1.

Figure 1 shows the coarsescale parameter obtained for the two reactions given in Table 1. Using this methodology, the coarsescale parameter as a function spatial mesh size can be obtained for a specific solidsolid reaction system. This parameter then can be used for largescale numerical simulation of a reactive system. Figure 1 shows that for both reactions the coarse parameter starts increasing about a dimensionless grid size of 5 suggesting that . Figure 2 shows the reaction rate as a function of dimensionless distance. The dimensionless size of the reaction zone is approximated by the region, where the reaction rate is larger than 0.01 times of the peak reaction rate. Results show that for both reactions the approximate dimensionless size of the reaction zone is about 13 suggesting that using finescale model one needs to use grid size of 5/13 of the size of the reaction zone to roughly capture the reaction front. The presented coarsescale formulation allows us to use larger grid size while maintaining accuracy of the solution. Figures 3 and 4 show comparisons of concentration and temperature profiles with and without use of the presented coarsescale formulation. Results show that the coarsescale formulation could represent the location of the reaction front accurately with a small number of grid cells. Figures 3 and 4 also show that by using coarsescale formulation one might choose a coarse spatial mesh size two times of the reaction zone while maintaining the numerical solution accuracy.
The calculated numerical errors for the reactions given in Table 1 are presented in Figure 5 for concentration and temperature for different grid sizes. Results show that the presented coarsescale formulation could significantly reduce the numerical error. In addition, the ratio of CPU time for coarsescale formulation and fine grid reference solutions is presented in Figure 6. Results show that the CPU time ratio scales with the inverse of grid size, suggesting a tenfold reduction in CPU time with a tenfold increase in grid size. Currently, we are working toward finding scaling or proportionality constant for multidimensional problems. CPU time savings is expected to be much more significant for multidimensional problems.
(a)
(b)
(c)
(d)
5. Concluding Remarks
Propagation of solidsolid reaction fronts often results in a thin reaction zone that is difficult to resolve numerically unless a large number of numerical grid cells are used. Such numerical simulations are computationally expensive to perform. In this study, a coarsescale formulation for numerical modeling of a onedimensional solidsolid reacting system is presented. The presented formulation is based on a Taylor series expansion of the reaction term and presents the modification of the reaction term in the coarse model that would allow an accurate solution. A key parameter in this formulation is or socalled coarsescale parameter, which is a function of the coarsescale spatial mesh size. Using the presented formulation, this parameter as a function of spatial mesh size can be obtained. This parameter then can be used for numerically solving a large scale and computationally intensive reacting system. It is shown that this formulation could reasonably obtain the accuracy of a fine grid numerical solution. It is shown that the coarsescale formulation considerably reduces the numerical error. In addition, it is revealed that the ratio of CPU times of a coarsescale model to that of a fine grid solution scales with the inverse of grid size. Such inverse proportionality implies a significant reduction in CPU time of a coarsescale model compared to that of a fine gridscale model for largescale simulations. Results obtained in this study for a onedimensional reacting system are promising in terms of reducing computational time. However, the presented methodology has a number of limitations that are the subject of our current research. First, the coarsescale parameter needs to be characterized as a function of the governing dimensionless numbers. In addition, the applicability of this method needs to be tested for multidimensional problems, where we expect a major reduction in CPU time. Presently, we are working toward implementation of the presented methodology for multidimensional problems and determination of the coarsescale parameter as a function of dimensionless numbers for solidsolid reactions. We anticipate that the computational time saving for multidimensional problems is more promising. The presented formulation improves our capabilities for conducting more accurate and faster numerical simulation of industrialscale solidsolid reactions.
Nomenclature
:  Heat capacity, 
:  Concentration, kg/m^{3} 
:  Molecular diffusion coefficient, m^{2}/s 
:  Activation energy, 
:  Heat of reaction, J/kg 
:  Preexponential rate constant, 
:  Gas constant, 
:  Time, s 
:  Temperature, K 
:  Coordinate, m. 
:  Coarsescale parameter 
:  Inverse of dimensionless activation energy 
:  Inverse of dimensionless heat of reaction 
:  Numerical error 
:  Coarsescale volume, m^{3} 
:  Finescale variable can be temperature, concentration, or reaction rate 
:  Reaction rate, 
:  Effective thermal conductivity, 
:  Density, 
:  Dimensionless temperature. 
a:  Adiabatic 
:  Critical 
:  Dimensionless 
:  Ignition 
0:  Initial value 
*:  Scale value. 
':  Deviation from coarse scale 
:  Average or coarse scale. 
Acknowledgments
The authors would like to acknowledge constructive comments from Dr. Brian Wood and Dr. Christof Meile. Helpful discussion with Dr. Mohsen Sadeghi is also acknowledged. The financial support of the Alberta Ingenuity Centre for In Situ Energy (AICISE) is acknowledged. The authors would like to thank the reviewers for useful comments.
References
 F. M. Lea, The Chemistry of Cement and Concrete, Edward Arnolds, London, UK, 3rd edition, 1970.
 N. M. Chechilo, R. J. Khvilivitskii, and N. S. Enikolopyan, “On the phenomenon of polymerization reaction spreading,” Doklady Akademii Nauk SSSR, vol. 204, pp. 1180–1181, 1972. View at: Google Scholar
 V. M. Ilyashenko, S. E. Solovyov, and J. A. Pojman, “Theoretical aspects of selfpropagating reaction fronts in condensed medium,” AIChE Journal, vol. 41, no. 12, pp. 2631–2636, 1995. View at: Publisher Site  Google Scholar
 W. D. Kingery, Introduction to Ceramics, John Wiley & Sons, New York, NY, USA, 1967.
 A. C. A. M. Bleijenberg, B. C. Lippens, and G. C. A. Schuit, “Catalytic oxidation of 1butene over bismuth molybdate catalysts. I. The system ${\text{Bi}}_{2}{\text{O}}_{3}{\text{MoO}}_{3}$,” Journal of Catalysis, vol. 4, no. 5, pp. 581–585, 1965. View at: Publisher Site  Google Scholar
 S. S. Tamhankar and L. K. Doraiswamy, “Analysis of solidsolid reactions: a review,” AIChE Journal, vol. 25, no. 4, pp. 561–582, 1979. View at: Publisher Site  Google Scholar
 Ya. B. Zeldovich and D. A. FrankKamenetzkii, “A theory of thermal propagation of flame,” ACTA PhysicoChimica URSS, vol. 9, no. 2, pp. 341–350, 1938. View at: Google Scholar
 S. B. Margolis, “An asymptotic theory of condensed twophase flame propagation,” SIAM Journal on Applied Mathematics, vol. 43, no. 2, pp. 351–369, 1983. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. B. Margolis, “An asymptotic theory of heterogeneous condensed combustion,” Combustion Science and Technology, vol. 43, no. 34, pp. 197–215, 1985. View at: Publisher Site  Google Scholar
 S. B. Margolis and R. C. Armstrong, “Two asymptotic models for solid propellant combustion,” Combustion Science and Technology, vol. 47, no. 12, pp. 1–38, 1986. View at: Publisher Site  Google Scholar
 J. Puszynski, J. Degreve, and V. Hlavacek, “Modeling of exothermic solidsolid noncatalytic reactions,” Industrial and Engineering Chemistry Research, vol. 26, no. 7, pp. 1424–1434, 1987. View at: Publisher Site  Google Scholar
 J. Rajaiah, H. Dandekar, J. Puszynski, J. Degreve, and V. Hlavacek, “Study of gassolid, heterogeneous, exothermic, noncatalytic reactions in a flow regime,” Industrial and Engineering Chemistry Research, vol. 27, no. 3, pp. 513–518, 1988. View at: Publisher Site  Google Scholar
 H. Dandekar, J. A. Puszynski, J. Degreve, and V. Hlavacek, “Reaction front propagation characteristics in noncatalytic exothermic gassolid systems,” Chemical Engineering Communications, vol. 92, pp. 199–224, 1990. View at: Publisher Site  Google Scholar
 B. D. Wood and S. Whitaker, “Diffusion and reaction in biofilms,” Chemical Engineering Science, vol. 53, no. 3, pp. 397–425, 1998. View at: Publisher Site  Google Scholar
 B. D. Wood and S. Whitaker, “Erratum to “Diffusion and reaction in biofilms”,” Chemical Engineering Science, vol. 55, no. 12, p. 2349, 2000. View at: Publisher Site  Google Scholar
 B. D. Wood and S. Whitaker, “Multispecies diffusion and reaction in biofilms and cellular media,” Chemical Engineering Science, vol. 55, no. 17, pp. 3397–3418, 2000. View at: Publisher Site  Google Scholar
 B. D. Wood, K. Radakovich, and F. Golfier, “Effective reaction at a fluidsolid interface: applications to biotransformation in porous media,” Advances in Water Resources, vol. 30, no. 67, pp. 1630–1647, 2007. View at: Publisher Site  Google Scholar
 C. Meile and K. Tuncay, “Scale dependence of reaction rates in porous media,” Advances in Water Resources, vol. 29, no. 1, pp. 62–71, 2006. View at: Publisher Site  Google Scholar
 A. G. Merzhanov, A. K. Filonenko, and I. P. Borovinskaya, “New phenomena in combustion of condensed systems,” Doklady Akademii Nauk SSSR, vol. 208, pp. 892–894, 1973. View at: Google Scholar
 A. G. Merzhanov and I. P. Borovinskaya, “A new class of combustion processes,” Combustion Science and Technology, vol. 10, no. 56, pp. 195–201, 1975. View at: Publisher Site  Google Scholar
 A. G. Merzhanov and B. I. Khaikin, “Theory of combustion waves in homogeneous media,” Progress in Energy and Combustion Science, vol. 14, no. 1, pp. 1–98, 1988. View at: Publisher Site  Google Scholar
 B. D. Wood, F. Cherblanc, M. Quintard, and S. Whitaker, “Volume averaging for determining the effective dispersion tensor: closure using periodic unit cells and comparison with ensemble averaging,” Water Resources Research, vol. 39, no. 8, p. 1210, 2003. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2009 Hassan Hassanzadeh et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.