Stiffness, numerical methods

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

Ordinaiy differential Eqs. (13-149) to (13-151) for rates of change of hquid-phase mole fractious are uouhuear because the coefficients of Xi j change with time. Therefore, numerical methods of integration with respect to time must be enmloyed. Furthermore, the equations may be difficult to integrate rapidly and accurately because they may constitute a so-called stiff system as considered by Gear Numerical Initial Value Problems in Ordinaiy Differential Equations, Prentice Hall, Englewood Cliffs, N.J., 1971). The choice of time... 

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. 

Existence and uniqueness of the particular solution of (5.1) for an initial value y° can be shown under very mild assumptions. For example, it is sufficient to assume that the function f is differentiable and its derivatives are bounded. Except for a few simple equations, however, the general solution cannot be obtained by analytical methods and we must seek numerical alternatives. Starting with the known point (tD,y°), all numerical methods generate a sequence (tj y1), (t2,y2),. .., (t. y1), approximating the points of the particular solution through (tQ,y°). The choice of the method is large and we shall be content to outline a few popular types. One of them will deal with stiff differential equations that are very difficult to solve by classical methods. Related topics we discuss are sensitivity analysis and quasi steady state approximation. 

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 fundamental difficulty in solving DEs explicitly via finite formulas is tied to the fact that antiderivatives are known for only very few functions / M —> M. One can always differentiate (via the product, quotient, or chain rule) an explicitly given function f(x) quite easily, but finding an antiderivative function F with F x) = f(x) is impossible for all except very few functions /. Numerical approximations of antiderivatives can, however, be found in the form of a table of values (rather than a functional expression) numerically by a multitude of integration methods such as collected in the ode... m file suite inside MATLAB. Some of these numerical methods have been used for several centuries, while the algorithms for stiff DEs are just a few decades old. These codes are... 

As an example, if only quasi-steady flow elements are used with volume pressure elements, a model s smallest volume size (for equal flows) will define the timescale of interest. Thus, if the modeler inserts a volume pressure element that has a timescale of one second, the modeler is implying that events which happen on this timescale are important. A set of differential equations and their solution are considered stiff or rigid when the final approach to the steady-state solution is rapid, compared to the entire transient period. In part, numerical aspects of the model will determine this, but also the size of the perturbation will have a significant impact on the stiffness of the problem. It is well known that implicit numerical methods are better suited towards solving a stiff problem. (Note, however, that The Mathwork s software for real-time hardware applications, Real-Time Workshop , requires an explicit method presumably in order to better guarantee consistent solution times.)... 

Numerical Method. Both the isothermal FFB and the adiabatic MAT models are very stiff due to the coke deactivation terms >j. The spline orthogonal collocation technique was used to solve the above models (19). Typically, the distance x was divided into two regions (0collocation points in each region. At the interface between the two regions, both the concentration and the mass flux were taken as continuous. The value of y varied with the degree of stiffness. 

The preceding concept of stability is not sufficient when stiff problems are considered and it is necessary to introduce the concept of absolute stability. A numerical method is said to be absolutely stable if the global discretization error remains bounded for a given step size h when the number, N, of steps tends to infinity. 

S. H. Lam,Singular Perturbations for Stiff Equations using Numerical Methods, in Recent Advances in Aerospace Sciences, C. Casci, ed.. New York Plenum Publishing Corp., 1985, 3-19. 

For the numerical study of the whole spectrum (for g R fixed), [79] uses a spectral tau-Chebychev discretization in y and the Arnold method (see [88]) to solve the generalized eigenvalue problem (see [89]). This numerical method is based on the orthonormalization of the Krylov space of the iterates of the inverse of the matrix A B. This method has been used more recently in [90]. It has been proven efficient in the stiff problems arising in the study of spectral stability of viscoelastic fluids. 

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. 

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. 

In section 2.2.4, stop condition was used to predict the maximum yield in a chemical reaction. In section 2.2.5, a stiff problem was solved using Maple s default numerical solver. We concluded that the conventional numerical methods fail for this stiff problem. Maple s stiff solver was found to be superior for this stiff problem. Generally, one has to use a stiff solver only if the conventional methods fail, as stiff solvers take more time to solve ordinary IVPs than the conventional solvers. 

Numerical Method of Lines for Stiff Nonlinear PDEs... 

In section 5.2.4, a stiff nonlinear PDE was solved using numerical method of lines. This stiff problem was handled by calling Maple s stiff solver. The temperature explodes after a certain time. The numerical method of lines (NMOL) technique was then extended to coupled nonlinear parabolic PDEs in section 5.2.5. By comparing with the analytical solution, we observed that NMOL predicts the behavior accurately. 

Both analytical and numerical methods of lines are presented in this chapter for elliptic partial differential equations. Semianalytical method, presented in this chapter is very powerful technique, and is valid for elliptic Partial differential equations with mixed boundaries also (Subramanian and White, 1999). Numerical method of lines presented in this chapter should be used with precaution, as it may not work for stiff problems. A total of seven examples were presented in this chapter. 

There are also other numerical methods for solving differential equations, which we do not discuss. The numerical methods can be extended to sets of simultaneous differential equations such as occur in the analysis of chemical reaction mechanisms. Many of these sets of equations have a property called stiffness that makes them difficult to treat numerically. Techniques have been devised to handle this problem, which is beyond the scope of this book. ... 

Thus, the composition of such a chemical system may be obtained, at time t, by solving the set of non-linear differential equations (76). The solution is generally obtained by numerical integration, which must be highly stable ones, as the problems encountered are locally exponential , giving rise to the phenomenon called. .stiffness . Warner gives a critical review of the different numerical methods... 

Numerical methods used to solve a system of ODEs are widely available in computational libraries and through texts such as Numerical Recipes. Certain considerations arise in the use of these standard techniques for nonlinear systems, particularly in models of chemical systems, which often consist of systems of stiff equations that require special care. Stiff equations are characterized by the presence of widely differing time scales, which leads to eigenvalues of the Jacobian matrix differing by many orders of magnitude. 

The equations in Table 2.6 can be solved using several numerical methods for ODEs. The resulting system is commonly stiff since the moments of the living polymer and active site concentration vary much faster than the other dependent variables in the system. A common way to eliminate this nuisance is to make the steady-state approximation for the moments of living polymer and active site concentration. This possibility will be examined in the following subsection. 

Cash, J. (2003). Efficient numerical methods for the solution of stiff initial-value problems and differential algebraic equations. Proceedings of the Royal Society, 459, 797-815. 

To determine the load distribution in a bolted connection accurately, taking into account the elastic properties of the joined members and the fasteners explicitly, it would be necessary to use a numerical method such as the finite element method (FEM). To greatly simplily the calculation of fastener load distribution it is assumed that the members are macroscopically rigid and that the elasticity is limited to local areas in the vicinity of the fasteners. It may furthermore be assumed that the load versus deformation response of an individual fastener is linear. Thus, the effect of member stiffness on fastener load distribution is taken into account by means of special correction factors. 

A major limitation of the simple method of determining the fastener load distribution is that the stiffness (elasticity and geometry) of the joined members is not explicitly taken into account. Another limitation of the simple method is that the specific elasticity properties of the analysed member is not considered when determining the fastener hole stress distribution. These limitations can effectively be overcome by adopting a numerical method such as the finite element method (FEM) for evaluation of the load distribution. 

Find solution to modelling equations using appropriate numerical methods this provides property characteristics of the unit cell at the current structural level. Estimate the effective stiffness matrix using one of the homogenization methods (Klusemann and Svendsen, 2010). 

---