Stiff Model Problems

The simplest linear model problem to represent the broad ideas in chemical kinetics is 

This equation, sometimes called the test equation in texts on numerical differential equations [13], has an important resemblance to chemical kinetics. Specifically, the rate of disappearance of y is proportional to y itself. As X (i.e., the rate constant) increases, the shorter the characteristic reaction time. The general solution to this problem is obviously 

The solution always reaches a steady state of y = 0, with X determining how fast it gets there. 

Regardless of the value of y, the characteristic time scale of this model equation is r = /X. Even at long time, when the solution is not changing at all (i.e., y = 0), the equation itself still has a characteristic time scale that can be quite short if X is large. Stiffness occurs in regions where the solution is changing slowly (or not at all), yet the characteristic time scales are very small. 

Consider the following first-order linear model problem [215], which has been designed to illustrate some of the computational issues associated with stiffness  

---

As a bit of an aside, one can think of the algebraic constraint as an infinitely stiff problem. Referring to the stiff model problem (Section 15.2), stiff problems are characterized by a fast transient and a slowly varying solution. Regardless of the initial condition, a stiff problem will always decay to the slowly varying solution, and the stiffer the problem, the faster will be the decay (e.g., Fig. 15.1). The situation in a problem like that in Fig. 7.5 is that there is no transient in the y2 component because it is a constraint, and not a differential equation. If, however, the y2 equation is modeled as y 2 = — X(y2 — 1), then y2 = (y2(0) — l)e Xl. As A. becomes larger, the differential equation becomes stiffer, and as X —> oo, the differential equation becomes an algebraic constraint. 

As our first model problem, we take the motion of a diatomic molecule under an external force field. For simplicity, it is assumed that (i) the motion is pla nar, (ii) the two atoms have equal mass m = 1, and (iii) the chemical bond is modeled by a stiff harmonic spring with equilibrium length ro = 1. Denoting the positions of the two atoms hy e 71, i = 1,2, the corresponding Hamiltonian function is of type... 

The standard discretization for the equations (9) in molecular dynamics is the (explicit) Verlet method. Stability considerations imply that the Verlet method must be applied with a step-size restriction k < e = j2jK,. Various methods have been suggested to avoid this step-size barrier. The most popular is to replace the stiff spring by a holonomic constraint, as in (4). For our first model problem, this leads to the equations d... 

In the model problem described earlier, replace the stiff spring potential Vs = K r - 1)2/2 by... 

The combination of time marching and Newton s method can be illustrated via a very simple model problem [277]. Consider two reactions, R + A B + P and R + B 2P, where in the first a reactant R reacts with a compound A to produce a compound B and a product P. Then, in the second reaction, R further reacts with B to produce two moles of P. If the reaction rates are significantly different, this will lead to a stiff system. For the sake of our example, presume that the mole fraction of R is fixed at a value of 0.1, and that the rate constants for the reactions are k = 1011 and ki = 1012, respectively. Furthermore take the equilibrium constants for the two reactions to be AT] =5 and K.2 = 15. With these parameters set, the mole fractions of A and B (A and B) are governed by the following system of equations. (The value of P is determined from the fact that the mole fractions must sum to unity.)... 

So-called stiff differential equation models are particularly challenging to solve. Stiff models have dynamic behavior that encompasses a wide range of time scales. An example would be fast kinetics combined with long fluid-residence times in a chemical reactor. Gear s method is perhaps the most commonly used technique for solving these types of problems. 

If none of the aforementioned information is available, then construction and inspection of models can give a rough idea of the angle strain involved. For example, cyclodecyne can be built with stiff models without problems, but cyclo-nonyne and cyclooctyne can only be constructed with a flexible C—C=C—C moiety. 

The third issue to be introduced is on the online model class selection. The usual approach in system identificatimi is to find the best/optimal model in a prescribed class of models, e.g., class of shear building models with uncertain interstory stiffnesses. This problem is commonly referred to as parametric identification. The more general problem of model class selection has not been explored as intensively as parametric identification. It is well known that a more complicated model class often fits the data better than one which has fewer adjustable parameters. 

The principle of the perfectly-mixed stirred tank has been discussed previously in Sec. 1.2.2, and this provides essential building block for modelling applications. In this section, the concept is applied to tank type reactor systems and stagewise mass transfer applications, such that the resulting model equations often appear in the form of linked sets of first-order difference differential equations. Solution by digital simulation works well for small problems, in which the number of equations are relatively small and where the problem is not compounded by stiffness or by the need for iterative procedures. For these reasons, the dynamic modelling of the continuous distillation columns in this section is intended only as a demonstration of method, rather than as a realistic attempt at solution. For the solution of complex distillation problems, the reader is referred to commercial dynamic simulation packages. 

The modeling of complex solids has greatly advanced since the advent, around 1960, of the finite element method [196], Here the material is divided into a number of subdomains, termed elements, with associated nodes. The elements are considered to consist of materials, the constitutive equations of which are well known, and, upon change of the system, the nodes suffer nodal displacements and concomitant generalized nodal forces. The method involves construction of a global stiffness matrix that comprises the contributions from all elements, the relevant boundary conditions and body and thermal forces a typical problem is then to compute the nodal displacements (i. e., the local strains) by solving the system K u = F, where K is the stiffness matrix, u the... 

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. 

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

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

Restrictions which may exist for the choice of a commercial reactor need not be imposed at the development stage. In some cases, a reactor of one type may be best for acquiring data in model characterisation, whereas a reactor of another type might be more suitable for full-scale production. (The cautions expressed in Sect. 4 must be taken into account.) Continuous flow back-mixed reactors can be very useful for kinetic studies because the absence of concentration gradients can reduce uncertainties in concentration measurements. When these reactors have attained a steady state, many of the problems associated with stiffness (see above) can be avoided. 

The first two sections of Chapter 5 give a practical introduction to dynamic models and their numerical solution. In addition to some classical methods, an efficient procedure is presented for solving systems of stiff differential equations frequently encountered in chemistry and biology. Sensitivity analysis of dynamic models and their reduction based on quasy-steady-state approximation are discussed. The second central problem of this chapter is estimating parameters in ordinary differential equations. An efficient short-cut method designed specifically for PC s is presented and applied to parameter estimation, numerical deconvolution and input determination. Application examples concern enzyme kinetics and pharmacokinetic compartmental modelling. 

Stiffness occures in a problem if there are two or more very different time scales on which the dependent variables are changing. Since at least one component of the solution is "fast", a small step size must be selected. There is, however, also a "slow" variable, and the time interval of interest is large, requiring to perform a large number of small steps. Such models are common in many areas, e.g., in chemical reaction kinetics, and solving stiff equations is a challenging problem of scientific computing. 

---