Improving Accuracy of Coarse Grid Numerical Solution of Solid-Solid Reactions by Taylor Series Expansion of the Reaction Term

Abstract

Exothermic solid-solid reactions lead to sharp reaction fronts that cannot be captured by coarse spatial mesh size numerical simulations that are often required for large-scale simulations. We present a coarse-scale 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 fine-scale behavior of a constant pattern reaction front while using a smaller number of numerical grid cells. Results for a one-dimensional solid-solid reacting system reveal reasonable computational time saving. The presented formulation improves our capabilities for conducting fast and accurate numerical simulations of industrial-scale solid-solid 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 solid-solid reactions. Modeling of solid-solid reactions is of significant interest in various industrial operations, including oxidation of metallic and nonmetallic mixed powders, cement industry [1], ferrites manufacturing, solid-state 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 reaction-transport from pore scale to continuum scale [7–18].

Accurate numerical simulation of solid-solid reactions is a challenging task due to the multiscale nature of the physical phenomena. Physical processes involved in solid-solid 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 subdiffusive-scale concentration and temperature gradients, while heat and mass diffusive processes have scales orders of magnitude larger than the reactions. Large-scale 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 industrial-scale simulation of the involved processes with high accuracy. In this paper, we provide a coarse-scale formulation to capture the subgrid-scale phenomena appropriate for large-scale numerical simulations. The paper is organized as follows. First, the fine-scale mathematical model used in this study is presented. Next, the coarse-scale model is described. Then, application of the coarse-scale model is given for a constant pattern reaction front, followed by summary and conclusions.

2. Fine-Scale Conservation Equations

The exothermic solid-solid reaction is assumed to take place in a one-dimensional semi-infinite 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 𝑡=0, the initial temperature 𝑇₀ at 𝑥=0 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 𝜕𝜃𝜕𝑡_𝐷=𝜕²𝜃𝜕𝑥²_𝐷+𝐶_𝐷𝜃exp(𝛽𝜃+1),𝜕𝐶_𝐷𝜕𝑡_𝐷=−𝛾𝐶_𝐷𝜃exp(𝛽𝜃+1), (2.1) where the following dimensionless groups are used [11, 13]: 𝜃=𝐸(𝑇−𝑇^∗)𝑅𝑇^∗2,𝐶_𝐷=𝐶𝐶₀,𝑡_𝐷=𝑡𝑡^∗,𝑡^∗=𝜌𝑐_𝑝𝑅𝑇^∗2𝐸𝑘(−Δ𝐻)𝐶₀𝐸exp(𝑅𝑇^∗𝑥),_𝐷=𝑥𝑥^∗,𝑥^∗=𝜆𝑡^∗𝜌𝑐_𝑝,𝛾=𝜌𝑐_𝑝𝑅𝑇^∗2𝐸(−Δ𝐻)𝐶₀,𝛽=𝑅𝑇^∗𝐸, (2) where T is temperature, C is concentration of reactant, 𝜌 is density, 𝑐_𝑝 is heat capacity, 𝜆 is the average thermal conductivity, k is the pre-exponential 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 𝜃=𝜃₀,0≤𝑥_𝐷≤∞,𝑡_𝐷𝐶=0,_𝐷=1,0≤𝑥_𝐷≤∞,𝑡_𝐷=0,𝜃=0,at𝑥_𝐷=0,𝑡_𝐷>0,𝜕𝜃𝜕𝑥_𝐷=0,𝑥→∞,𝑡_𝐷>0,𝜕𝐶_𝐷𝜕𝑥_𝐷=0,𝑥→∞,𝑡_𝐷>0. (2.3) 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. Coarse-Scale 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 coarse-scale variable can be defined by 1𝜓=𝑣_𝑣𝜓𝑑𝑣, (3.1) 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 coarse-scale numerical model can be expressed by [18] 𝜕𝜃𝜕𝑡_𝐷=𝜕²𝜃𝜕𝑥²_𝐷+𝜕ℜ,𝐶_𝐷𝜕𝑡_𝐷=−𝛾ℜ, (3.2) where ℜ=𝐶_𝐷𝜃exp(𝛽𝜃+1),1ℜ=𝑣_𝑣𝐶_𝐷𝜃exp(𝛽𝜃+1)𝑑𝑣=𝐶_𝐷𝜃exp()𝛽𝜃+1 (3.3) are the fine-scale reaction rate and coarse-scale average reaction rate, respectively. The reaction rate ℜ can be approximated by a Taylor series expansion around the coarse-scale temperature and concentration. The above expansion can be represented in terms of deviation from the coarse-scale variables as given by ℜ≈ℜ(𝐶_𝐷,𝜃)+𝜕ℜ𝜕𝐶_𝐷|_{𝐶𝐷,𝜃}𝐶^′_𝐷+𝜕ℜ|𝜕𝜃_{𝐶𝐷,𝜃}𝜃^′+1(𝜕2!²ℜ𝜕𝐶²_𝐷|_{𝐶𝐷,𝜃}(𝐶^′_𝐷)²+2𝜕²ℜ𝜕𝐶_𝐷|𝜕𝜃_{𝐶𝐷,𝜃}𝐶^′_𝐷𝜃^′+𝜕²ℜ𝜕𝜃²|_{𝐶𝐷,𝜃}(𝜃^′)²)+⋯, (3.4) 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 coarse-scale reaction rate can be expressed by ℜ≈ℜ(𝐶_𝐷,1𝜃)+(𝜕2!²ℜ𝜕𝐶²_𝐷|_{𝐶𝐷,𝜃}(𝐶^′_𝐷)²+2𝜕²ℜ𝜕𝐶_𝐷|𝜕𝜃_{𝐶𝐷,𝜃}𝐶^′_𝐷𝜃^′+𝜕²ℜ𝜕𝜃²|_{𝐶𝐷,𝜃}(𝜃^′)²𝜕)+⋯,₂ℜ𝜕𝜃²|_{𝐶𝐷,𝜃}=1−2𝛽(1+𝛽𝜃)(𝛽𝜃+1)⁴𝐶_𝐷exp(𝜃𝛽𝜕𝜃+1),₂ℜ𝜕𝐶²_𝐷|_{𝐶𝐷,𝜃}𝜕=0,₂ℜ𝜕𝐶_𝐷|𝜕𝜃_{𝐶𝐷,𝜃}=exp(𝜃/(𝛽𝜃+1))(𝛽𝜃+1)², (3.5)(3.6)(3.7)(3.8) where by definition the terms with first derivatives are dropped in the averaging process [18]. Equation (3.5) is a coarse-scale representation of the reaction rate that includes both coarse-scale 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 coarse-scale concentration and temperature. Instead, we intend to express the quantities (𝜃^′)² and 𝐶^′_𝐷𝜃^′ in a form proportional to 𝜕𝜃/𝜕𝑥_𝐷 and (𝜕𝜃/𝜕𝑥_𝐷)(𝜕𝐶_𝐷/𝜕𝑥_𝐷) such that the final coarse-scale model can be presented as a function of coarse-scale temperature and concentration only. Similar approach has been used by Meile and Tuncay [18] for linear convection-diffusion-reaction 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 fine-scale and the coarse-scale 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 (𝜃^′)²𝜕=𝛼(𝜃𝜕𝑥_𝐷)²,𝐶^′_𝐷𝜃^′𝜕=𝛼(𝜃𝜕𝑥_𝐷𝜕)(𝐶_𝐷𝜕𝑥_𝐷). (3.9) By substituting (3.9) in (3.5), the coarse-scale equations (3.2) can be written as 𝜕𝜃𝜕𝑡_𝐷=𝜕²𝜃𝜕𝑥²_𝐷+ℜ(𝐶_𝐷,𝛼𝜃)+(2!2𝜕²ℜ𝜕𝐶_𝐷|𝜕𝜃_{𝐶𝐷,𝜃}(𝜕𝜃𝜕𝑥_𝐷𝜕)(𝐶_𝐷𝜕𝑥_𝐷𝜕)+²ℜ𝜕𝜃²|_{𝐶𝐷,𝜃}(𝜕𝜃𝜕𝑥_𝐷)²𝜕),𝐶_𝐷𝜕𝑡_𝐷=−𝛾ℜ(𝐶_𝐷,𝛼𝜃)−𝛾{(2!2𝜕²ℜ𝜕𝐶_𝐷|𝜕𝜃_{𝐶𝐷,𝜃}(𝜕𝜃𝜕𝑥_𝐷𝜕)(𝐶_𝐷𝜕𝑥_𝐷𝜕)+²ℜ𝜕𝜃²|_{𝐶𝐷,𝜃}(𝜕𝜃𝜕𝑥_𝐷)²)}, (3.10) where 𝛼 is the proportionality constant or coarse-scale parameter and is a function of coarse grid size. The coarse-scale 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 right-hand 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 coarse-scale formulation can be rendered as 𝜕𝜃𝜕𝑡_𝐷=𝜕²𝜃𝜕𝑥²_𝐷+ℜ(𝐶_𝐷,𝛼𝜃)+𝜕2!²ℜ𝜕𝜃²|_{𝐶𝐷,𝜃}(𝜕𝜃𝜕𝑥_𝐷)²,𝜕𝐶_𝐷𝜕𝑡_𝐷=−𝛾ℜ(𝐶_𝐷,𝛼𝜃)−𝛾{+𝜕2!²ℜ𝜕𝜃²|_{𝐶𝐷,𝜃}(𝜕𝜃𝜕𝑥_𝐷)²}. (3.11) 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 ℜ(𝐶_𝐷,𝜃)+𝛼(𝜕²ℜ/𝜕𝜃²)_{𝐶𝐷,𝜃}(𝜕𝜃/𝜕𝑥_𝐷)²/2.

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 coarse-scale model predictions from the fine-scale reference solution, we define a numerical error using the following expression as a measure of accuracy: 𝜀={[∫(𝜓−𝜓)²𝑑𝑥_𝐷][∫𝜓²𝑑𝑥_𝐷]}^1/2, (3.12) 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 coarse-scale formulation (3.11) is complete and can be used for subsequent large-scale simulation of a solid-solid reactive system. In the following, we use the above procedure for determining the parameter 𝛼 as a function of coarse spatial mesh size for some solid-solid reacting systems.

4. Application of the Coarse-Scale Formulation

Based on the scaling groups used in nondimensionalizing the governing equations, the frontal behavior of a solid-solid reaction depends on two parameters, namely, 𝛽 (inverse of dimensionless activation energy) and 𝛾 (inverse of dimensionless heat of reaction). The coarse-scale formulation described previously is applied for two solid-solid reactions. The governing differential equations are discretized using an explicit-in-time finite difference approximation. A block-centered 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 large-scale 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 solid-solid 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 (1/𝛾). Figure 1 shows 𝛼 as a function of numerical spatial mesh size for different solid-solid 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]: 𝛾_𝑐=12(𝜃_𝑐−𝜃₀), (4.1) where 𝜃_𝑐=𝜃₀+21−2𝛽. (4.2) For 𝛾<𝛾_𝑐, a reacting system demonstrates constant pattern reaction front propagation. The parameter 𝛾_𝑐 for reactions studied is given in Table 1.

Table 1 

Dimensionless parameters for solid-solid reacting systems used in this study.

Figure 1 

Coarse-scale parameter 𝛼 as a function of grid size for two reactions in Table 1.

Figure 1 shows the coarse-scale parameter 𝛼obtained for the two reactions given in Table 1. Using this methodology, the coarse-scale parameter 𝛼as a function spatial mesh size can be obtained for a specific solid-solid reaction system. This parameter then can be used for large-scale 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 𝑥^∗=0.2Δ𝑥. 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 fine-scale 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 coarse-scale 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 coarse-scale formulation. Results show that the coarse-scale 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 coarse-scale formulation one might choose a coarse spatial mesh size two times of the reaction zone while maintaining the numerical solution accuracy.

Figure 2 

Reaction rate for the two reactions as a function of dimensionless distance.

Figure 3 

Concentration and temperature distributions at 𝑡𝐷=400 for Reaction 1 given in Table 1 with and without coarse-scale formulation for three grid sizes of 25, 16.7, and 9.1. (—) reference solution, () without using coarse-scale parameter, and () with using coarse-scale parameter.

Figure 4 

Concentration and temperature distributions at 𝑡𝐷=650 for Reaction 2 given in Table 1 with and without coarse-scale formulation for three grid sizes of 25, 16.7, and 9.1. (—) reference solution, () without using coarse-scale parameter, and () with using coarse-scale parameter.

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 coarse-scale formulation could significantly reduce the numerical error. In addition, the ratio of CPU time for coarse-scale 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 ten-fold reduction in CPU time with a ten-fold 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)

Figure 5 

Numerical error in temperature (left) and concentration (right) distributions for reactions given in Table 1 with and without coarse-scale formulation for different grid sizes.

Figure 6 

Ratio of coarse to fine grid (reference) CPU times versus dimensionless grid size for the two reactions given in Table 1.

5. Concluding Remarks

Propagation of solid-solid 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 coarse-scale formulation for numerical modeling of a one-dimensional solid-solid 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 so-called coarse-scale parameter, which is a function of the coarse-scale 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 coarse-scale formulation considerably reduces the numerical error. In addition, it is revealed that the ratio of CPU times of a coarse-scale 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 coarse-scale model compared to that of a fine grid-scale model for large-scale simulations. Results obtained in this study for a one-dimensional 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 coarse-scale 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 coarse-scale parameter as a function of dimensionless numbers for solid-solid 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 industrial-scale solid-solid reactions.

Nomenclature

𝑐𝑝:	Heat capacity, Jkg−1K−1
𝐶:	Concentration, kg/m3
𝐷:	Molecular diffusion coefficient, m2/s
𝐸:	Activation energy, Jkmol−1
Δ𝐻:	Heat of reaction, J/kg
𝑘:	Pre-exponential rate constant, s−1
𝑅:	Gas constant, Jkmol−1K−1
𝑡:	Time, s
𝑇:	Temperature, K
𝑥:	Coordinate, m.

Greek Letters

𝛼:	Coarse-scale parameter
𝛽:	Inverse of dimensionless activation energy
𝛾:	Inverse of dimensionless heat of reaction
𝜀:	Numerical error
𝑣:	Coarse-scale volume, m3
𝜓:	Fine-scale variable can be temperature, concentration, or reaction rate
ℜ:	Reaction rate, kgm−3s−1
𝜆:	Effective thermal conductivity, Jm−1s−1K−1
𝜌:	Density, kgm−3
𝜃:	Dimensionless temperature.

Subscripts

a:	Adiabatic
𝑐:	Critical
𝐷:	Dimensionless
𝑖:	Ignition
0:	Initial value
*:	Scale value.

Superscripts

':	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.