Stretching-based diagnostics and reduction of chemical kinetic models with diffusion

Abstract

A new method for diagnostics and reduction of dynamical systems and chemical kinetic models is proposed. The method makes use of the local structure of the normal stretching rates by projecting the dynamics onto the local directions of maximal stretching. The approach is computationally very simple as it implies the spectral analysis of a symmetric matrix. Notwithstanding its simplicity, stretching-based analysis derives from a geometric basis grounded on the pointwise applications of concepts of normal hyperbolicity theory. As a byproduct, a simple reduction method is derived, equivalent to a “local embedding algorithm”, which is based on the local projection of the dynamics onto the “most unstable and/or slow modes” compared to the time scale dictated by the local tangential dynamics. This method provides excellent results in the analysis and reduction of dynamical systems displaying relaxation towards an equilibrium point, limit cycles and chaotic attractors. Several numerical examples deriving from typical models of reaction/diffusion kinetics exhibiting complex dynamics are thoroughly addressed. The application to typical combustion models is also analyzed.

Introduction

Many chemical, biochemical and biological processes involve a large number of species and reactions [1], [2], and the resulting kinetic schemes, either under perfectly mixed condition or, a fortiori when spatial inhomogeneities (diffusion and/or convection) are accounted for, are expressed by means of high-dimensional dynamical systems. The occurrence of a large number of different time scales, and the resulting stiffness of the model equations are common features of all these models.

In order to highlight the role and the influence of the different kinetic steps (model diagnostics), simplify the model equations, and possibly obtain a reduced model, different theoretical and computational approaches have been proposed, such as the Intrinsic Low-Dimensional Manifold (ILDM) method [3], [4], [5], Computational Singular Perturbation (CSP) [6], [7], the Method of Invariant Manifold (MIM) [8], [9], the method by Fraser and Roussel [10], [11], methods based on Lyapunov functions such as the thermodynamic free energy [12], [13], methods based on the intrinsic dynamics in the tangent/cotangent bundle [14], the Zero Time Derivative method [15], a combination of the Fraser and Roussel method with CSP [16]. For a review and a comparison of several of these approaches see [17], [18].

The common denominator of all these computational strategies is a simple and evocative paradigm: the occurrence within the phase space of a Slow Invariant Manifold (SIM) for system dynamics, attracting nearby orbits, possessing a lower dimensionality than the phase space, and representing the backbone around which orbit dynamics is organized.

In the case of infinite-dimensional systems, such as reaction–diffusion models [19] expressed by means of a system of partial differential equations, the transposition of ILDM and CSP approaches have been proposed by Hadjinicolau and Goussis [20], Singh et al. [21], Goussis et al. [22], via the concept of an infinite-dimensional slow manifold obtained by projecting the dynamics onto the local slow modes associated with the kinetics. These approaches have been developed aside from the inertial manifold theory, proposed by Temam et al. [23], [24] that provides functional–theoretical criteria and computational methods [25], [26], [27], [28] to obtain a finite-dimensional representation of reaction–diffusion dynamics within an absorbing invariant set.

Motivated by the case study of singularly perturbed systems [29], [30], [31], considered as a benchmarking prototype, the recent literature on model simplification and reduction has been characterized by a superposition of tools and methods deriving from geometric concepts (essentially expressing manifold invariance), perturbation theory and computational strategies. CSP by Lam and Goussis [6], [7], and the approach proposed in [16], [32] provide examples of the superposition of paradigms.

Although the connection between singularly perturbed system and computational methods for model reduction in generic models is strong and motivated, it is important to observe that the formulation of reduction strategies for generic dynamical systems of physical interest is in general divorced from a perturbative analysis, and that singularly perturbed systems represent solely a nonexhaustive class of dynamical models which can be tackled and simplified by means of model reduction methods.

The intermingling of concepts deriving from singular perturbation theory, differential geometry and numerical analysis has provided a wealth of different computational strategies to tackle model simplification of complex kinetics schemes. On the other hand, some controversial issues have been raised on the definition of slow manifolds and their properties. The reason for this may be attributed to different reasons: (i) the conditions imposed by different authors on slow manifolds may have different nature (e.g. by imposing some smoothness and analyticity criteria on the local representation of the manifold itself that by other authors are neglected [33], [34], [35]); (ii) any geometrical definition of the slow manifold should be grounded on global properties defined throughout the entire phase space, including the behavior at infinity [36], while this may not be the case in perturbation studies; (iii) in practical applications to model reduction of complex kinetic schemes, “intrinsic low-dimensional manifolds” may lack some basic properties (such as invariance [37], see e.g. [38]), and this limitation “collides” with more formal mathematical definitions of slow invariant manifolds [39]. Correspondingly, the very basic concept of “slow/fast decomposition” of complex reaction schemes may involve some intrinsic degree of arbitrariness since it relies on the specific method adopted in model diagnostics and reduction (such as ILDM [3], CSP [7], MIM [8]).

Apart from these fundamental issues on the definition of slow invariant manifolds, it is important to observe that most of the computational strategies proposed for model simplification and reduction have been developed and tested for chemically reacting systems evolving towards a stable equilibrium point. The ILDM and CSP methods have been verified and benchmarked by considering chemical reacting systems and combustion models relaxing towards equilibria as paradigmatic examples [4], [40]. It is, therefore, not surprising that the application of these model reduction techniques display problems and shortcomings when applied to dynamical systems evolving towards more general limit sets such as limit cycles and chaotic attractors [14]. In point of fact, the occurrence of persistent asymptotic oscillations defined on periodic/aperiodic/chaotic limit sets is a common feature observed in many chemical and biochemical open systems (closed adiabatic chemical systems obeying the law of mass action, possess a unique stable equilibrium point). For a review, see the monographs by Scott and Goldbeter [41], [42], and references cited therein.

As far as reduction strategies are concerned, the occurrence of more complex limit attractors introduces an extra degree of complexity to be taken into account. An efficient reduction algorithm for these dynamical systems should not only correctly describe the relaxation towards the limit set, but also the oscillating behavior on the invariant limit set itself. Several geometrical approaches have been proposed for tackling this wider class of dynamical systems, which are based on invariant vector dynamics within the tangent bundles [14], [43], [44]. All these approaches are theoretically extremely interesting, since they address invariant geometric features of dynamical systems on a fundamental level, but suffer the problem of being computationally onerous for higher-dimensional dynamical systems, as they require the explicit estimate of vector dynamics along system orbits. In the case one is interested in the reconstruction of the dynamics exclusively on invariant limit sets, a wealth of different approaches have been provided. For a review, see e.g. [45] and references cited therein.

The aim of this article is to propose a simple and efficient model reduction and diagnostic method for generic dynamical systems regardless of their finite or infinite-dimensional nature, and of the geometry of their limit sets. At the core of the method is a geometric characterization based on local normal stretching rates. The numerics requires the spectral characterization of a symmetric matrix.

In spite of its conceptual simplicity, the stretching-based method for system diagnostics is based on a geometric description of local tangent and normal dynamics. This geometric description finds its theoretical justification in the theory of normal hyperbolicity [39], [46], viewed on a local level, i.e., pointwisely along system orbits. To some extent, the stretching-based method shares some analogies with some recent modification of ILDM [8], [47]. However, there are some conceptual differences between the method proposed and these ILDM-based strategies, as discussed in the remainder of the article (Section 3).

The stretching-based reduction method can be viewed as a local embedding technique stemming from the stretching rate analysis of the system, obtained by locally projecting the dynamics onto the most unstable/slow directions. The reduction method proposed is computationally simple and efficient as it exclusively involves the solution of a lower-dimensional system of ordinary differential equations without nonlinear constraints, the dimension of which coincides with the number of relevant normal directions that should be locally accounted for.

The article is organized as follows. Section 2 reviews succinctly the basic notions of normal hyperbolicity and the definition of local normal/tangential stretching rates. Several examples taken from prototypical kinetic models highlight the meaning of the stretching-based analysis in the characterization of dynamic properties along system orbits and invariant manifolds. Section 3 introduces the stretching-based approach to model characterization in the general n-dimensional case and explores connections and differences with other existing methods. In particular, we show how to compute the normal stretching rate spectrum (and the corresponding set of directions of maximum normal stretching restricted to the normal subspace) and how it can be used for performing a local classification of the slow and fast (or unstable/stable) modes of the dynamics. Section 4 addresses the stretching-based method via some paradigmatic examples, by considering both low-dimensional models and nonlinear reaction–diffusion kinetics originating periodic and chaotic oscillations. The role of conservation laws in chemical kinetics is addressed in an Appendix, by considering the Michaelis–Menten enzymatic reaction as a prototypical model. Section 5 develops the stretching-based reduction method, and discusses several empirical criteria for defining the number of relevant normal modes. Apart from low-dimensional prototypical systems exhibiting chaotic behavior, the examples of application of the stretching-based reduction strategy concern infinite-dimensional reaction–diffusion systems and premixed flames.