Numerical Solution of Stiff Equations

The successful numerical solution of differential equations requires attention to two issues associated with error control through time step selection. One is accuracy and the other is stability. Accuracy requires a time step that is sufficiently small so that the numerical solution is close to the true solution. Numerical methods usually measure the accuracy in terms of the local truncation error, which depends on the details of the method and the time step. Stability requires a time step that is sufficiently small that numerical errors are damped, and not amplified. A problem is called stiff when the time step required to maintain stability is much smaller than that that would be required to deliver accuracy, if stability were not an issue. Generally speaking, implicit methods have very much better stability properties than explicit methods, and thus are much better suited to solving stiff problems. Since most chemical kinetic problems are stiff, implicit methods are usually the method of choice. 

The mathematical models of the reacting polydispersed particles usually have stiff ordinary differential equations. Stiffness arises from the effect of particle sizes on the thermal transients of the particles and from the strong temperature dependence of the reactions like combustion and devolatilization. The computation time for the numerical solution using commercially available stiff ODE solvers may take excessive time for some systems. A model that uses K discrete size cuts and N gas-solid reactions will have K(N + 1) differential equations. As an alternative to the numerical solution of these equations an iterative finite difference method was developed and tested on the pyrolysis model of polydispersed coal particles in a transport reactor. The resulting 160 differential equations were solved in less than 30 seconds on a CDC Cyber 73. This is compared to more than 10 hours on the same machine using a commercially available stiff solver which is based on Gear s method. 

For a system of S chemical species and R reactions c is the S vector of concentrations, k the R vector of time independent parameters (rate coefficients), and f the vector of the R rate expression functions. If the overall reaction is isothermal and takes place in a well-mixed vessel, equation (1) comprises a detailed chemical kinetic model (DCKM) of the reaction. The integration of the model equations can present difficulties because the rate coefficients may vary from one another by many orders of magnitude, and the differential equations are stiff. Numerical methods for the solution of stiff equations are discussed by Kee et al. [1]. Efficient solvers for stiff sets of equations have been developed and are available in various software packages. Some of these are described in Chapter 5. Additional information can be found in Refs. [2,3]. 

Steihaug, T. and Wolfbrandt, A., An Attempt to Avoid Exact Jacobian and Non-linear Equations in the Numerical Solution of Stiff Differential Equations, Math. Comp. 33, 521-534 (1979). 

The numerical solution of the equation was achieved by using a library routine which selected the order and step size automatically and was designed for stiff differential equations (IMSL-DGEAR). The routine was applied in the following iterative scheme ... 

Another way of classifying the integration techniques depends on whether or not the method is explicit, semi-implicit, or implicit. The implicit and semi-implicit methods play an important role in the numerical solution of stiff differential equations. To maintain the continuity of the section, we will first describe the explicit integration techniques in the context of one-step and multistep methods. The concept of stiffness and implicit methods are considered in a separate subsection, which also marks the end of this section. 

Families of finite elements and their corresponding shape functions, schemes for derivation of the elemental stiffness equations (i.e. the working equations) and updating of non-linear physical parameters in polymer processing flow simulations have been discussed in previous chapters. However, except for a brief explanation in the worked examples in Chapter 2, any detailed discussion of the numerical solution of the global set of algebraic equations has, so far, been avoided. We now turn our attention to this important topic. 

As mentioned earlier, overall accuracy of finite element computations is directly detennined by the accuracy of the method employed to obtain the numerical solution of the global system of algebraic equations. In practical simulations, therefore, computational errors which are liable to affect the solution of global stiffness equations should be carefully analysed. 

Check the difficulty of numerical solution due to equation stiffness. 

Dynamic simulations are also possible, and these require solving differential equations, sometimes with algebraic constraints. If some parts of the process change extremely quickly when there is a disturbance, that part of the process may be modeled in the steady state for the disturbance at any instant. Such situations are called stiff, and the methods for them are discussed in Numerical Solution of Ordinary Differential Equations as Initial-Value Problems. It must be realized, though, that a dynamic calculation can also be time-consuming and sometimes the allowable units are lumped-parameter models that are simplifications of the equations used for the steady-state analysis. Thus, as always, the assumptions need to be examined critically before accepting the computer results. 

This again shows that the component balance may thus have different time constants, which depend on the relative magnitudes of the equilibrium constants Ki which again can lead to problems of numerical solution due to equation stiffness. 

A commercial stiff ordinary differential equation solver subroutine, DVOGER, is available in the IMSL Library (3). This subroutine uses Gear s method for the solution of stiff ODE s with analytic or numerical Jacobians. The pyrolysis model was solved using DVOGER and the analytical Jacobians of Eqs. (14) and (15). For a residence time of 0.0511 in dimensionless time, defined as t/t where 9... 

Warner [176] has given a comprehensive discussion of the principal approaches to the solution of stiff differential equations, including a hundred references among the most pertinent books, papers and application packages directed at simulating kinetic models. Emphasis has been put not only on numerical and software problems such as robustness, improving the linear equation solvers, using sparse matrix techniques, etc., but also on the availability of a chemical compiler, i.e. a powerful interface between kineticist and computer. 

Over the past ten years the numerical simulation of the behavior of complex reaction systems has become a fairly routine procedure, and has been widely used in many areas of chemistry, [l] The most intensive application has been in environmental, atmospheric, and combustion science, where mechanisms often consisting of several hundred reactions are involved. Both deterministic (numerical solution of mass-action differential equations) and stochastic (Monte-Carlo) methods have been used. The former approach is by far the most popular, having been made possible by the development of efficient algorithms for the solution of the "stiff" ODE problem. Edelson has briefly reviewed these developments in a symposium volume which includes several papers on the mathematical techniques and their application. [2]... 

Q = [a +iaRe —c)Y/ . For boundary layer instability problems, i e —> 00 and then Q >> laj. This is the source of stiffness that makes obtaining the numerical solution of (2.3.21) a daunting task. This causes the fundamental solutions of the Orr-Sommerfeld equation to vary by different orders of magnitude near and far away from the wall. This type of behaviour makes the governing equation a stiff differential equation that suffers from the growth of parasitic error, while numerically solving it. 

Given the initial conditions (concentrations of the 22 chemical species at t = 0), the concentrations of the chemical species with time are found by numerically solving the set of the 22 stiff ordinary differential equations (ODE). An ordinary differential equation system solver, EPISODE (17) is used. The method chosen for the numerical solution of the system includes variable step size, variable-order backward differentiation, and a chord or semistationary Newton method with an internally computed finite difference approximation to the Jacobian equation. 

One problem with the availability of a large set of standard libraries is that at times it is possible to use the wrong module for a given task. A common situation involves the solution of differential equations. Some systems of differential equations contain derivatives that vary over wide scales, and these are known as stiff systems of differential equations. Therefore, a stiff differential equation solver should be used in these cases otherwise, substantial numerical errors or convergence problems will result. 

---