Stiff ODE Models

With this transformation, the normal equations now become, 

As we have already pointed out in this chapter for systems described by algebraic equations, the introduction of the reduced sensitivity coefficients results in a reduction of cond(A). Therefore, the use of the reduced sensitivity coefficients should also be beneficial to non-stiff systems. 

We strongly suggest the use of the reduced sensitivity whenever we are dealing with differential equation models. Even if the system of differential equations is non-stiff at the optimum (when k=k ), when the parameters are far from their optimal values, the equations may become stiff temporarily for a few iterations of the Gauss-Newton method. Furthermore, since this transformation also results in better conditioning of the normal equations, we propose its use at all times. This transformation has been implemented in the program for ODE systems provided with this book. 

---

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

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. 

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. 

Depending on the numerical techniques available for integration of the model equations, model reformulation or simplified version of the original model has always been the first step. In the Sixties and Seventies simplified models as sets of ordinary differential equations (ODEs) were developed. Explicit Euler method or explicit Runge-Kutta method (Huckaba and Danly, 1960 Domenech and Enjalbert, 1981 Coward, 1967 Robinson, 1969, 1970 etc) were used to integrate such model equations. The ODE models ignored column holdup and therefore non-stiff integration techniques were suitable for those models. 

However, in batch distillation, the system is frequently very stiff, owing either to wide ranges in relative volatilities or large differences in tray and reboiler holdups. Therefore, if methods for non-stiff problems are applied to stiff problems (ODE models but having column holdup and/or energy balances), a very small integration step must be used to ensure that the solution remains stable (Meadow, 1963 Distefano, 1968 Boston et al., 1980 Holland and Liapis 1983, etc.). 

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

Dabdub, D., and Seinfeld, J. H. (1995) Extrapolation techniques used in the solution of stiff ODEs associated with chemical kinetics of air quality models, Atmos. Environ. 29, 403-410. 

The model for the schematic system (Fig. 3) consists of the simple ODEs (3) (or (6)), (7), (8) and (9), which form an initial value problem (IVP). In the case that pure hydrogen is used, its pressure is kept constant and the liquid-phase components are nonvolatile, the gas-phase balance equations (8)-(9) are discarded and the gas-phase concentration in eqs (3) and (6) is obtained e.g. from the ideal gas law. The initial conditions, i.e. the concentrations at t=0 are equal everywhere in the system and the IVP can be solved numerically by any stiff ODE-solver. 

The governing parabolic PDFs describing the model were discretized with respect to the spatial coordinates of the catalyst particles and the column length coordinate. The resulting ODEs were solved numerically with the sparse version of the stiff ODE solver, LSODES (Appendix 2). 

Often, models are considered to be simple , if it is an easy task to set up their equations of motion, or if there are few equations which perhaps are nearly linear. Being simple in that sense does not necessarily imply that the integration task is easy. We will see this when discussing stiff ODEs. 

The partial-differential equations (3) are solved numerically through finite difference approximation for the spatial derivatives and the method of line for time advancement. The model medium is represented by a discretized line with a resolution from 50 up to 200 points. The resulting set of ordinary differential equations is integrated with a stiff ODE solver [85]. Care is taken to vary the spatio-temporal resolution in order to check the reliability of the reported phenomena. 

In summary, in the equilibrium-chemistry limit, the computational problem associated with turbulent reacting flows is greatly simplified by employing the presumed mixture-fraction PDF method. Indeed, because the chemical source term usually leads to a stiff system of ODEs (see (5.151)) that are solved off-line, the equilibrium-chemistry limit significantly reduces the computational load needed to solve a turbulent-reacting-flow problem. In a CFD code, a second-order transport model for inert scalars such as those discussed in Chapter 3 is utilized to find ( ) and and the equifibrium com-... 

The Runge-Kutta algorithm cannot handle so-called stiff problems. Computation times are astronomical and thus the algorithm is useless, for that class of ordinary differential equations, specialised stiff solvers have been developed. In our context, a system of ODEs sometimes becomes stiff if it comprises very fast and also very slow steps and/or very high and very low concentrations. As a typical example we model an oscillating reaction in The Belousov-Zhabotinsky (BZ) Reaction (p.95). 

You may wonder why we would ever be satisfied with anything less than a very accurate integration. The ODEs that make up the mathematical models of most practical chemical engineering systems usually represent a mixture of fast dynamics and slow dynamics. For example, in a distillation column the liquid flow or hydraulic dynamic response occurs fairly rapidly, of the order of a few seconds per tray. The composition dynamics, the rate of change of hquid mole fractions on the trays, are usually much slower—minutes or even hours for columns with many trays. Systems with this mixture of fast and slow ODEs are called stiff systems. 

Despite its success, the embedded model approach still requires repeated solution of the process model (and sensitivities). For large processes or for processes that require the solution of rigorous underlying procedures, this approach can become expensive. Moreover, for stiff or otherwise difficult systems, this approach is only as reliable as the ODE solver. The embedded model approach also offers only indirect ways of handling time-dependent constraints. Finally, the optimal solution of this approach is only as good as its control variable parameterization, which often can only be improved by a priori information about the specific problem. Consequently, we now consider the simultaneous approach to (16) as an alternative to solution methods for (17). 

The kinetic profiles displayed in Figure 7.16 have been integrated numerically with MATLAB s stiff solver ode 15s using the rate constants k, = 1000 M s k2 = 100 s 1 for the initial concentrations [A]0 = 1M, Cut ) = 10 4M, and /1 0 = [C]0 = 0 M. For this model, the standard Runge-Kutta routine is far too slow and thus useless. In MATLAB Example 7.7b, the function is given that generates the differential equations. It is repeatedly called by the ODE-solver in MATLAB Example 7.7a. 

---