Stiff ODEs

Stiffness often refers to a system with a very short time constant. For the system model 

Even if a is very large, if s t) is a slowly varying function, the system responds slowly. To follow this slow response numerically is problematic since the stability of the method depends only on a and has nothing to do with s(t). A very small time step might be required, even though the overall system response is a slowly varying function. This is one example of a phenomenon called stiffness. 

Another example of stiffness occurs with a system (two or more) of ODEs each with greatly different time constants. For example, consider the system 

The very short time constant in the second equation requires a very small h to follow both y and z in time. Problems of instability can easily arise. To solve stiff problems, the most typical approach is to use an implicit method, which is known to exhibit excellent stability properties (recall the backward Euler method). Special software packages are available for solving stiff systems. Fortunately, many of the straightforward chemical engineering problems enconntered in practice do not yield stiff systems, but when difficulties arise, stiffness might well be the culprit. 

The standard Euler methods and Runge-Kutta methods do not converge for stiff ODE S. A still system can be defined as one in which the stability of the numerical methods used becomes an issue. Maple has an inbuilt stiff solver. 

---

When the dynamic system is described by a set of stiff ODEs and observations during the fast transients are not available, generation of artificial data by interpolation at times close to the origin may be very risky. If however, we ob-... 

ODE solver. Relative to non-stiff ODE solvers, stiff ODE solvers typically use implicit methods, which require the numerical inversion of an Ns x Ns Jacobian matrix, and thus are considerably more expensive. In a transported PDF simulation lasting T time units, the composition variables must be updated /Vsm, = T/At 106 times for each notional particle. Since the number of notional particles will be of the order of A p 106, the total number of times that (6.245) must be solved during a transported PDF simulation can be as high as A p x A sim 1012. Thus, the computational cost associated with treating the chemical source term becomes the critical issue when dealing with detailed chemistry. 

Fractional time stepping is widely used in reacting-flow simulations (Boris and Oran 2000) in order to isolate terms in the transport equations so that they can be treated with the most efficient numerical methods. For non-premixed reactions, the fractional-time-stepping approach will yield acceptable accuracy if A t r . Note that since the exact solution to the mixing step is known (see (6.248)), the stiff ODE solver is only needed for (6.249), which, because it can be solved independently for each notional particle, is uncoupled. This fact can be exploited to treat the chemical source term efficiently using chemical lookup tables. 

In a transported PDF simulation, the chemical source term, (6.249), is integrated over and over again with each new set of initial conditions. For fixed inlet flow conditions, it is often the case that, for most of the time, the initial conditions that occur in a particular simulation occupy only a small sub-volume of composition space. This is especially true with fast chemical kinetics, where many of the reactions attain a quasi-steady state within the small time step At. Since solving the stiff ODE system is computationally expensive, this observation suggests that it would be more efficient first to solve the chemical source term for a set of representative initial conditions in composition space,156 and then to store the results in a pre-computed chemical lookup table. This operation can be described mathematically by a non-linear reaction map ... 

However, in practice, the integral must be evaluated using a stiff ODE solver or chemical lookup tables (see Section 6.9). Because transported PDF simulations are typically used for reacting flows with complex chemistry, the chemical-reaction step will often dominate the overall computational cost. It is thus important to consider carefully the computational efficiency of the chemical-reaction step when implementing a transported PDF simulation. 

The situation is different for those readers who do not have access to Matlab and rely completely on Excel. In the following, we explain how a fourth order Runge-Kutta method can be incorporated into a spreadsheet and used to solve non-stiff ODE s. 

P. Seifert, Computational experiments with algorithms for stiff ODEs, Computing, 38 (1987) 163-176. 

In what follows, we begin by introducing two examples of process systems with recycle and purge. First, we analyze the case of a reactor with gas effluent connected via a gas recycle stream to a condenser, and a purge stream used to remove the light impurity present in the feed. In the second case, the products of a liquid-phase reactor are separated by a distillation column. The bottoms of the column are recycled to the reactor, and the trace heavy impurity present in the feed stream is removed via a liquid purge stream. We show that, in both cases, the dynamics of the system is modeled by a system of stiff ODEs that can, potentially, exhibit a two-time-scale behavior. 

An Iterative Approach for the Solution of a Class of Stiff ODE Models of Reacting Polydispersed Particles... 

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. 

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

These results show that the proposed technique provides a fast and reliable method for the solution of stiff ODE models of reacting polydispersed particles. Recently, Turton (9) applied this method successfully to the modeling of wood char combustion in a transport reactor. 

One particular class of ordinary differential equation solvers (ODE-solvers) handles stiff ODEs and these are widely known as stiff solvers. In our context, a system of ODEs sometimes becomes stiff if it comprises very fast and also very slow steps or relatively high and low concentrations. A typical example would be an oscillating reaction. Here, a highly sophisticated step-size control is required to achieve a reasonable compromise between accuracy and computation time. It is well outside the scope of this chapter to expand on the intricacies of modem numerical... 

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]... 

Implementation of the calculation followed in general the scheme of Dougherty et al. [8] The kinetic problem Eq. (1) was solved separately using the BELLCHEM code consisting of a chemical compiler followed by a stiff ODE solver using the method of Gear. [//] Results were stored on disc for subsequent input to the sensitivity calculation. 

The numerical Method of Lines as implemented in the routine NDSolve of the Mathematica system deals with system (32) by employing the default fourth order finite difference discretization in the spatial variable Z, and creating a much larger coupled system of ordinnary equations for the transformed dimensionless temperature evaluated on the knots of the created mesh. This resulting system is internally solved (still inside NDSolve routine) with Gear s method for stiff ODE systems. Once numerical results have been obtained and automatically interpolated by NDSolve, one can apply the inverse expression (31.b) to obtain the full dimensionless temperature field. 

A possible way to achieve this is to solve equation (4.1), together with all the sensitivity equations belonging to different parameters. In this case, the stiff ode solver has to decompose a large (m + l)n x (m + l)n Jacobian at each time-step, which is very inefficient. If the kinetic ode (4.1) is coupled with a single sensitivity equation at one time, the joint Jacobian is smaller, but the kinetic system of odes has to be solved m times. 

The most efficient algorithm for the solution of the sensitivity differential equations is called the decoupled direct method (ddm), which was first applied in chemical kinetics by Dunker [67, 68]. He drew attention to the fact that equations (4.1) and (4.6) have the same Jacobian, so that a stiff ode solver will use the same step size and order of approximation in the solution of both odes. The ddm method first takes a step for the solution of equation (4.1) and then performs steps for the solution of equation (4.6) for / = 1,. . . , m. The procedure is repeated in the subsequent steps. Since the Jacobian of the equations is the same, it has to be triangularized only once for each time interval. This method is applied in the program SENKIN [69]. 

The most recent version of POLYMATH has both a normal and a stiff ODE solver along with a library to store home problems worked using POLYMATH. The example problems in the text that use POLYMATH are in the POLYMATH library in the CD-ROM. 

Simultaneous to the graph creation, kinetic properties in each vRxn are used to create the appropriate reaction rate equations (ordinary differential equations, ODE). These properties include rate constants (e.g., Michaelis constant, Km, and maximum velocity, Vmax, for enzyme-catalyzed reactions, and k for nonenzymatic reactions), inhibitor constants, A) and modes of inhibition or allosterism. The total set of rate equations and specified initial conditions forms an initial value problem that is solved by a stiff ODE equation solver for the concentrations of all species as a function of time. The constituent transforms for the each virtual enzyme are compiled by carefully culling the literature for data on enzymes known to act on the chemicals and chemical metabolites of interest. 

Using the boundary conditions (equations (5.54) and (5.55)) the boundary values uo and un+i can be eliminated. Hence, the method of lines technique reduces the nonlinear parabolic PDE (equation (5.48)) to a nonlinear system of N coupled first order ODEs (equation (5.52)). This nonlinear system of ODEs is integrated numerically in time using Maple s numerical ODE solver (Runge-Kutta, Gear, and Rosenbrock for stiff ODEs see chapter 2.2.5). The procedure for using Maple to solve nonlinear parabolic partial differential equations with linear boundary conditions can be summarized as follows ... 

Byrne, G.D., Hindmarsh, A.C. Stiff ODE Solvers A Review of Current and Coming Attractions. Journal of Computational Physics 70(1), 1-62 (1987)... 

---