Digital simulation of thermal reactions
Introduction
There is continuing interest in the mathematics describing thermal explosions, and solutions for various geometries. The physical system is often a gas confined within some bounds, reacting chemically so that it heats up and either attains a steady state or explodes. The steady state was described by Frank-Kamenetskii [1] who later also considered the time development of such a reaction [2], and already in 1939 gathered several parameters in the equation in what was to be named the Frank-Kamenetskii parameter , see below. At values below its critical value , a steady state is reached for a given geometry and boundary conditions, and an explosion ensues for values above it. The critical value is obviously of considerable interest.
The governing equation for the process is [2] (but using slightly different notation)where T is the temperature, the time, k the thermal diffusivity , c is the specific heat capacity , Q the heat of the reaction , z the pre-exponential Arrhenius parameter , E the activation energy of the reaction, R the gas constant, and X the spatial variable, expressing the distance from the centre of the body for either a slab with thickness 2a and infinite length, a cylinder with radius a and infinite length, or a sphere of radius a. Thus D specifies the dimension of the body in which the reaction takes place. It is commonly assumed that no substance is consumed by the reaction. Boundary conditions areThe boundary condition at was justified by Frank-Kamenetskii on the assumption that the heat capacity of the reacting substance (often a gas) is much smaller than that of the vessel containing it [2], and that assumption will be followed here. Other conditions at the boundary at are possible; Several authors [4], [6], [5], [3] considered Robin conditions at this boundary, which can be formulated asOur assumed boundary condition will then correspond to , so that (3) is reduced to (2b).
There have also been works that assume reacting substance being consumed during the reaction [7], [8], [3], but this will not be followed in this work.
It is convenient to normalise the variables:whereexpresses the rate of the reaction at ambient temperature , and is often set to zero; it will often be small (but its effect will be examined below). is the so-called Frank-Kamenetskii parameter [1], [9], [2]. There is a critical value of this parameter, , for a given geometry, below which the system attains a steady state and above which the system explodes (possibly, see below). The above normalisations lead to the dimensionless form of (1):with the new boundary conditionsThe present authors are not aware of analytical solutions to the above equation with its boundary conditions; that is, the time development of u. Frank-Kamenetskii [9], [2] derived a solution for the steady state for in the slab system ():in which a is the solution ofFor the cylindrical system a solution was derived by Rice [10] (for small x), Boddington and Gray [11], Harley [12], citing Frank-Kamenetskii [2] and Chambré [13]. Harley’s expression iswithHarley also provides an alternative (mathematically equivalent) solution. The present authors are not aware of an analytical solution for the spherical geometry.
Most studies have focussed on the existence of a steady state, described by (6) with the left-hand side set to zero, but the development of the temperature with time is also of interest. This has been studied in several works [14], [12], [8], [3]. In such a case, as there is no analytical solution, numerical solutions must be found. Harley [12] simulated the problem, developing time-dependent profiles . He used the hopscotch method [15], comparing it to what he calls the Crank–Nicolson method [16]; he concluded that both are roughly equal in efficiency, while (not surprisingly) the built-in pde solver of MATHEMATICA was much slower. However, Harley did not make use of the full potential of Crank–Nicolson, which essentially makes use of approximations to values halfway between time nodes by taking averages of the point values at the present time and the next. Doing this with (6) results in a nonlinear set of equations, which Harley avoided by discretising the chemical term explicitly. This however degrades the method and we will show below that Crank–Nicolson outperforms Harley’s implementation when the nonlinearity is dealt with in a style consistent with the Crank–Nicolson method.
Hopscotch was shown to have a significant disadvantage. It is a stable method, so that one, in principle, has a free choice of the parameter ,in which is the length of the time interval in the simulation and h the spatial interval. However, Feldberg pointed out [17] that hopscotch becomes increasingly inaccurate as is increased because, for large , changes near the boundaries do not propagate sufficiently into the interior of the diffusion space. For this reason, electrochemical simulators, who at one time favoured hopscotch for its ease of use [18], [19], [20], [21], have largely stopped using the method. The Crank–Nicolson method does have one drawback: it can produce large error oscillations, often damped only slowly, if there is an initial sharp transient in the field, such as the response to a step function. These can be damped [23], [22] but this is not needed here, as there is no sharp initial transient, so Crank–Nicolson is a useful and efficient method in the present context.
Harley also computed steady state temperatures, by prolonging the time marching simulation until no further changes were seen. These then matched the analytical values closely. It is however unnecessary to do this; the steady state can be computed directly by solving (1) or (6), after setting the left-hand side to zero. This is the method used in this work.
Section snippets
Steady state solutions – the effect of
For there are solutions only for , the critical Frank-Kamenetskii parameter, and in this range there exist two solutions for the steady state for the slab and cylinder, and multiple solutions for a range of in the spherical case. Physically it is the lower of the solutions that is attained (although Joseph [24] speculates on a physical upper solution state, albeit unstable, under certain conditions). For there exist solutions for arbitrary values of . However there is a critical
Discretisation of the time-marching pde (6)
Eq. (6) can be discretised in the Crank–Nicolson sense. We divide the range into N equal intervals of length h. At a given time level t we have values and wish to compute the values for , denoted as . (6) then becomesfor . The boundary condition (7b) gives .
The boundary condition (7c) at can be handled in two different ways. One can discretise it
Steady state computations
The steady state computations amount to solving a two-point boundary value problem for a second order ordinary differential equation. We do this by the method of ‘shooting’, i.e. we guess an initial value such that together with we have an initial value problem. At we use the above mentioned Maclaurin expansion such that the discretization of (18) becomesWith given and (21) can readily be solved for . Then every subsequent can be explicitly
Computational
Computations were performed at three separate locations and using separate programs in order to decrease the danger of error; (i) under Linux on an IBM x345 computer using FreePascal; (ii) under Windows 7 (64 bits) on a PC equipped with an AMD Athlon® II X4 processor, using the gfortran 95/2003 compiler V.4.6 and (iii) under Linux running on an Amitech Blueline PC and Intel Fortran 90/95 V.11.3. All computations were run using IEEE 754 Standard double precision, giving approximately 16-decimal
Critical
Table 1 shows the results of our computations of for a range of values. These are reliable to . The values for agree with those by others [14], [4], [11], [13], [1], [9], [2], [5], [33], [34], [35], most of whom present only a few decimals, but the value of Moise and Pritchard [34] for the sphere and , 3.32199212, agrees with ours to eight significant figures, and those of Boddington and Gray [11] for all three geometries (and ), given to five decimal places, also agree
Conclusions
Numerical solutions of the equation for thermal expansion of a reacting gas have been carried out, exploring both steady state and time-marching solutions. The critical Frank-Kamenetskii parameter has been evaluated to seven decimal places for the slab, cylinder and spherical geometry and the role of the critical activation parameter was explored. It was found that there exist one or more mathematical steady states for any if , the curves for steady temperature at the center of the
References (42)
- et al.
The dependence of criticality on activation-energy when reactant consumption is neglected
Combust. Flame
(1978) - et al.
Direct evaluation of branching points for equations arising in the theory of explosions of solid explosives
J. Comp. Phys.
(1975) - et al.
The numerical solution for reaction-diffusion combustion with fuel consumption
Appl. Math. Comp.
(2005) Hopscotch method: the numerical solution of the Frank-Kamenetskii partial differential equation
Appl. Math. Comp.
(2010)Propagational inadequacy of the hopscotch finite difference algorithm: the enhancement of performance when used with an exponentially expanding grid for simulation of electrochemical diffusion problems
J. Electroanal. Chem.
(1987)- et al.
Hopscotch: An algorithm for the numerical solution of electrochemical problems
J. Electroanal. Chem.
(1984) - et al.
Chronoamperometry at an ensemble of microdisk electrodes
J. Electroanal. Chem.
(1984) - et al.
Influence of insulation geometry on the current at microdisk electrodes
J. Electroanal. Chem.
(1984) - et al.
Explicit hopscotch and implicit finite-difference algorithms for the Cottrell problem: exact analytical results
J. Electroanal. Chem.
(1986) - et al.
Damping of Crank–Nicolson error oscillations
Comput. Biol. Chem.
(2003)
Non-linear heat generation and stability of the temperature distribution in conducting solids
Int. J. Heat Mass Trans.
Reactive-diffusive equation with variable pre-exponential factor
Mech. Res. Commun.
Calculation of thermal explosion limits
Acta Physiochim. URSS
Diffusion and Heat Transfer in Chemical Kinetics
Reaction-diffusion model for combustion with fuel consumption: Ii. Robin boundary conditions
IMA J. Appl. Math.
The method of fundamental solutions for non-linear thermal explosions
Commun. Numer. Meth. Eng.
Reaction-diffusion model for combustion with fuel consumption: I. Dirichlet boundary conditions
IMA J. Appl. Math.
Analytical solution of the problem of thermal explosions in a cylindrical enclosure
Zh. Fiz. Khim.
The role of heat conduction in thermal gaseous explosions
J. Chem. Phys.
Temperature profiles in endothermic and exothermic reactions and interpretation of experimental rate data
Proc. Roy. Soc. Lond Ser A - Mat. Phys. Sci.
On the solution of the Poisson–Boltzmann equation with application to the theory of thermal explosions
J. Chem. Phys.
Cited by (14)
Approximate analytical solution for mathematical models of thermal ignition and non-isothermal catalytic zero order reaction in a spherical geometry
2019, Journal of King Saud University - Engineering SciencesCitation Excerpt :Eq. (2) can also describe the steady state differential heat balance for catalyst particles in which non-isothermal zero order reaction takes place. Hlavacek and Marek (1968) used the Frank-Kamenetskii approximation (Eq. (3)) and presented the analytical solutions for the cases of solid and hollow slab and cylinder and the numerical solution for the case of sphere (Britz et al., 2011). Aris (1975) devoted a good part of his book to discuss the Lane-Emden equation and Frank-Kamenetskii equation and their relation to each other.
Approximate solution for the Lane-Emden equation of the second kind in a spherical annulus
2019, Journal of King Saud University - Engineering SciencesCitation Excerpt :Several theoretical studies and numerical methods were used for the study of the initial value problem of the isothermal Lane-Emden equation or the related boundary value problem of the Frank-Kamenetskii equation (Adler, 2011; Aris, 1975;Bazley and Wake, 1981,Britz et al., 2011; Chandrasekhar, 1967; Enig, 1967; Frank-Kamenetskii, 1955; Gustafson and Eaton, 1982; Hlavacek and Marek, 1968; Moise and Pritchard, 1989; Nazari-Golshan et al., 2013; Steggerda, 1965).
Flow, thermal criticality and transition of a reactive third-grade fluid in a pipe with Reynolds' model viscosity
2016, Journal of HydrodynamicsExothermic Diffusion-Reaction and Explosion Branch-Chain of an Asymmetrical Heating Channel with Convective Cooling
2023, Combustion Science and TechnologyOn Criticality for a Branched-chain Thermal Reactive-Diffusion in a Cylinder
2022, Combustion Science and Technology