The Rate-Controlled Constrained-Equilibrium Approach to Far-From-Local-Equilibrium Thermodynamics

Abstract

The Rate-Controlled Constrained-Equilibrium (RCCE) method for the description of the time-dependent behavior of dynamical systems in non-equilibrium states is a general, effective, physically based method for model order reduction that was originally developed in the framework of thermodynamics and chemical kinetics. A generalized mathematical formulation is presented here that allows including nonlinear constraints in non-local equilibrium systems characterized by the existence of a non-increasing Lyapunov functional under the system’s internal dynamics. The generalized formulation of RCCE enables to clarify the essentials of the method and the built-in general feature of thermodynamic consistency in the chemical kinetics context. In this paper, we work out the details of the method in a generalized mathematical-physics framework, but for definiteness we detail its well-known implementation in the traditional chemical kinetics framework. We detail proofs and spell out explicit functional dependences so as to bring out and clarify each underlying assumption of the method. In the standard context of chemical kinetics of ideal gas mixtures, we discuss the relations between the validity of the detailed balance condition off-equilibrium and the thermodynamic consistency of the method. We also discuss two examples of RCCE gas-phase combustion calculations to emphasize the constraint-dependent performance of the RCCE method.

1. Introduction

In thermodynamics, gas dynamics, chemical kinetics, control theory, information theory, economics, ecology, biology, and other sciences, it is increasingly important to determine the non-equilibrium time-dependent behavior of the state of a system, for which the description of the equilibrium states is known and given by the solution of a constrained maximization problem. The most accurate way of doing this consists of specifying the state variables, developing a detailed kinetic model (DKM) of the non-equilibrium system, obtaining a full set of rate equations for all the state variables based on the postulated DKM, and integrating this set of equations to yield the time-dependent behavior of the non-equilibrium state. Since this approach requires, in principle, modeling all the detailed mechanisms, a potential concern is to obtain the kinetic rate data for all the detailed mechanisms. However, for example in the context of chemical kinetics, such data are quite uncertain and involve a lot of guess work. Moreover, when the DKM method must be coupled to transport equations and turbulence models or interface dynamics, the modeling and computational tasks often become formidable because of the intrinsic presence of a wide range of time and length scales, which may result in the well-known stiffness difficulties.

In the framework of chemical kinetics, such difficulties have motivated the development of numerous model order reduction techniques during the past two decades [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26]. Moreover, for micro- and nano-systems in which the small numbers of particles entail that discreteness and stochasticity play an important role on the time evolution of the species population, similar reduction methods have been proposed within the growing field of Stochastic Chemical Kinetics [27,28,29] and Stochastic Simulation Algorithms [30]. One common feature shared by the aforementioned techniques is that they all need a large DKM as the starting point. However, even though theoretically the most accurate description of the kinetics, independent of the context, comes from the corresponding DKM, a general concern in modeling complex chemical systems for which the underlying DKM involves a large number of reaction mechanisms stems from the fact that the rates of a vast majority of reactions are unknown or, at best, highly uncertain.

An alternative method which can help reduce the modeling and computational efforts is the Rate-Controlled Constrained-Equilibrium (RCCE) method, in which the state of the system at any time is that which maximizes an entropy function or minimizes an appropriate form of free energy function subject to a small number of constraints, namely, state functionals. At each instant in time, the state of the system, which is described in terms of internal variables such as the reaction coordinates (in the chemical kinetics context), is assumed to be at equilibrium with respect to all the internal dynamical mechanisms that cannot alter the values of a small number of rate-controlling constraining functionals. The mechanisms that instead do contribute to changing the values of the constraints control the time evolution through their impact on the rate of change of the values of the constraining functionals. For a closed homogeneous reacting system, the resulting time evolution (trajectory in state space) is a continuous shifting of the assumed constrained equilibrium state towards higher entropies until the highest entropy (stable equilibrium) state is reached. For an open heterogeneous system subject to stationary boundary conditions, transport phenomena will also contribute to the rate of change of the local values of the constraints, and the continuous shifting of the assumed constrained local equilibrium states will evolve until the (stable) steady state is reached. The choice in RCCE can also be viewed as the best pragmatic choice according to Bayesian Inference in the framework of Jaynes’ Maximum Entropy Principle (see, e.g., [31,32,33]).

Two familiar examples of RCCE approach are the assumptions of locally complete and locally frozen chemical equilibrium. The former is characterized by atomic elements as the only constraints to the composition, subject to which pressure, temperature and Gibbs free energy are defined. The latter instead treat the entire chemical composition as a constant of the motion. Even in the absence of chemical reactions, another familiar example of the RCCE approach is the assumption of local thermodynamic equilibrium employed in the basic equations of fluid dynamics to describe non-equilibrium states by means of local temperature, pressure, and chemical potentials. Obviously, this assumption breaks down in the presence of fast diffusion processes, e.g., across a shock wave.

An application of RCCE reduces the total number of equations (algebraic + differential) to the total number of constraints, which for systems including heavy hydrocarbons can be dramatically smaller than the number of species in a corresponding DKM. Nonetheless, the dynamic evolution of all “left-out” species can be determined based on the requirement of constrained equilibrium. This often provides a good first order estimate of how kinetically important a species is when RCCE is used for model development. RCCE calculations have provided accurate results in relatively complex hydrocarbon combustion problems [34,35]. Two typical examples are discussed at the end of the paper. The effectiveness of the RCCE method and the accuracy of its results depend on the identification of a proper set of rate-controlling constraining functionals, which requires a fundamental understanding of the detailed dynamics of the system. In other words, the deeper the physical insight into the system’s dynamics, the more effective and structurally consistent a set of constraining functionals can be identified. In each particular field of application, a specific expertise must be developed in order to identify the controlling mechanisms and to group them into classes which must each be characterized by one or more functionals that are to become the rate-controlling constraints.

The objective in RCCE is to identify and group the rate limiting mechanisms into classes each characterized by a functional that is time-invariant within the class and which at least at some stage during the time evolution is slowly varying with respect to the full dynamics. In general, due to time-scale diversity and time-varying competition between different mechanisms, not all the constraining functionals are necessarily slowly varying and rate-controlling during the entire dynamic evolution. By identifying which constraining functionals are rate controlling during a given stage (for example, the ignition delay stage in a combustion process) it is possible to identify the corresponding underlying class of rate limiting mechanisms and gain useful physical insight.

Furthermore, despite what it may imply, a DKM is a reduced model itself for the simple reason that it neither includes all possible species in the reactor nor all the complex reactions governing them. From this point of view, we note that RCCE could be to a large extent independent of the level of approximation with which the rates of the less important reaction schemes are usually estimated in a DKM. The reason is that RCCE assumes at each instant in time the complete chemical equilibrium composition (highest entropy) compatible with the instantaneous local values of the state variables and the assumed constraining functionals. In other words, an informed application of RCCE puts these reactions in equilibrium, and it is well known that the complete chemical equilibrium composition is independent of the underlying detailed kinetic scheme ([36], p. 578). It is true that many of such reactions do not impose rate limiting constraints on the system’s evolution and are, therefore, eliminated in most DKM-based model order reduction techniques. But eliminating a reaction means freezing it. Instead, in RCCE the majority of such reactions are put in (constrained) equilibrium rather than being frozen.

The rates of change of the constraining functionals are to be estimated in principle from the full dynamics, but in practice only the faster mechanisms within the corresponding class are sufficient. In general, these faster mechanisms are often the most important in the fully detailed scheme and, therefore, they are the ones for which the most reliable experimental data on reaction rates are most likely available. Moreover, the physical insight, which is gained when a new successful constraining functional is identified and understood in connection with the corresponding underlying class of rate limiting mechanisms, is likely to be useful also to further improve, simplify, and complete the DKM.

In reacting gas dynamics, the RCCE method has been applied successfully for the description and simplification of many specific non-equilibrium problems. Sugden and co-workers [37], Kaskan [38], Schott [39], and Franciscus and Lezberg [40] treated the hydrogen-oxygen combustion problem. Keck and Gillespie [41] (see also [42,43]) developed a general rate-controlled constrained-equilibrium method for reacting gas mixtures. Delichatsios and Keck [44], Sarofim and co-workers [45], and Morr and Heywood [46] applied the method to predict carbon monoxide and nitric oxide formation in hydrocarbon combustion. Galant and Appleton [47] further developed the method for open reacting systems. Levine [48] treated reactions in molecular beams. Procaccia and Ross [49] investigated the stability in highly non-equilibrium systems. It has also formed the basis for Close Parallel Inertial Manifold (CPIM) [50] and Inertial Constrained-Equilibrium Pre-Image Curve (ICE-PIC) [51]. The method has been later developed and applied by several others, see e.g., [34,35,50,52,53,54,55,56,57,58,59].

Ideas similar to RCCE have been applied in specific applications also in other fields. For example, Kerner [60] studied speciation in ecological systems. Dornbush [61] developed a theory of exchange rate overshooting due to a changing monetary policy based on the notion that asset market prices including exchange rates adjust to the new equilibrium in practically no time whereas market prices of goods are “sticky”, adjust in finite time and thus rate-control the relaxation of the system towards the new equilibrium state.

It should be noted [62] that the entropy maximization done by RCCE at each instant of time does not necessarily imply that the relaxation is along the path of steepest entropy ascent [63,64,65,66]. This is consistent with the theoretical proof by Ziegler [67] as well as a recent numerical study [17] that maximal entropy production rate is not a general feature of the standard equations of chemical kinetics.

The objective of this paper is to present a generalized formulation of the RCCE method valid also for nonlinear constraints, together with a rigorous step-by-step derivation of its implementation in the standard framework of chemical kinetics of gas mixtures with linear constraints, a discussion about its general thermodynamic consistency, and examples about the physical meaning of several basic constraints that have proved effective in simplifying combustion modeling. The chemical kinetics implementation is presented with strong emphasis on the thermodynamics aspects, careful discussion of the assumptions, explicit indication of all functional dependences, and discussion of the most important features of the RCCE method, but the generalized formulation allows its application also in diverse nonequilibrium contexts, provided the time evolution is characterized by the existence of a Lyapunov functional which cannot be increased by the system’s internal dynamics.

The paper is structured as follows. In Section 2, we present a generalized formulation of the mathematical structure of the problem of describing the time evolution of non-equilibrium states. In Section 3, we discuss the general ideas behind the constrained-equilibrium method, and we effectively extend the method to the case of nonlinear constraints. In Section 4, we present a step-by-step detailed implementation of the foregoing general formulation in the framework of nonequilibrium thermodynamics of a mixture (ideal and nonideal) subject to complex chemical kinetics. There, by giving detailed proofs and explicit dependences, we address the subtleties about the domain of validity of detailed balance and the general nonnegativity of the rate of entropy generation. Finally, in Section 5 we discuss via two examples the physical meaning of the basic constraints that have proved most effective in reducing combustion modeling. Section 6 draws our conclusions.

2. General Non-Equilibrium Problem

In this and the next sections, we discuss the formulation of non-equilibrium dynamics for a generic state description which captures the model structure of a large class of systems of interest in various areas of science and engineering.

Let us start by asserting that when building a mathematical model of a system (physical, biological, chemical, etc.) we explicitly or implicitly select a “level of description” which postulates the existence of a set of elementary building blocks which, for the purposes of the model, we consider as indivisible and unalienable. For example, a given physical object may be described at various levels, depending on the phenomena we set out to model. In particular, as we further discuss below, the choice of the level of description is related to the range of time scales which characterize the phenomena of interest. So, if we are not interested in chemical and nuclear reactions, i.e., we may assume that on the time scale of interest they are practically “frozen”, then we may assume that the elementary constituents are the atoms and the molecules, as we do when we study non-reacting transport phenomena. If we are interested in chemical reactions but not nuclear reactions, we may take the elements (or, more precisely, the atomic nuclei and the electrons) as the elementary constituents. If we are interested in fission and fusion nuclear reactions, we may take neutrons, protons, and electrons. Or at a deeper level, we may take quarks and so on. Again, to describe resonance interactions between electromagnetic radiation and the electrons in atoms and molecules, it is in many instances sufficient to consider as elementary constituents the photons and a limited number of electronic levels.

Once the system’s constituents have been selected, the model defines what variables identify the different possible states. Let us denote by $x=\left\{x_1,x_2,\dots x_{n_{\mathrm{x}}}\right\}$ the vector of $n_{\mathrm{x}}$ variables describing the state of the system under study [68]. Often, it is convenient to separate the vector x into a vector of internal state variables that are not directly accessible to external control and a vector of externally controllable parameters, i.e., $x=\left\{x^{\mathrm{int}},x^{\mathrm{ext}}\right\}$. For example, in a closed chemically reacting system, species concentrations are internal variables controlled by the internal kinetics, energy and volume may be considered external when their changes are prescribed by heat and work interactions through the system’s boundaries and compressing or expanding displacements of the boundaries, respectively.

Next, let us assume that the dynamics of the system is described by a rate equation (equation of motion) of the form

$$\dot{x}=D^{\mathrm{int}}\left(x\right)+D^{\mathrm{ext}}(x,t)$$

(1)

where $D^{\mathrm{int}}\left(x\right)$ represents the internal or endogenous dynamical mechanisms and $D^{\mathrm{ext}}(x,t)$ the external or exogenous mechanisms. For example, if the system is a control volume in a reacting fluid flow, the endogenous mechanisms are those responsible for the spontaneous relaxation towards a local equilibrium state whereas the exogenous mechanisms are those accounting for the transport of energy, momentum and matter across the boundaries of the control volume.

The equilibrium states of the system are those that would not be changing in time if the system were isolated, i.e., if $D^{\mathrm{ext}}(x,t)=0$. Thus, a state $x^{\mathrm{eq}}$ is equilibrium if

$$D^{\mathrm{int}}\left(x^{\mathrm{eq}}\right)=0$$

(2)

In general, Equation (2) admits more than one solution, namely, the system admits more than one equilibrium state. For example, the second law of thermodynamics requires that a physical system has at least one equilibrium state for each of the possible values of the energy [69].

The system’s definition may impose elemental or geometrical constraints implying that a set of functionals that have values that are either fixed or externally controlled, i.e.,

$$C_i^{\infty ,\mathrm{str}}\left(x\right)=b_i$$

(3)

for all the states x of the system, where $b_k$ are constants or externally specified functions of time. These functionals generally specify restrictions on the geometry of the state space, or internal interconnections between internal partitions, or elemental conservation intrinsic in the chosen level of description. The reason for the superscript ∞ will become apparent in the next section. For example, if the state variables $x_i$ represent probabilities then a geometrical constraint is $C_1^{\infty ,\mathrm{str}}\left(x\right)={\sum }_i{x}_i$ with $b_1=1$. Again, if state variables $x_1=V_1$, $x_2=V_2$, $x_3=V_3$ represent the volumes of three regions of space in which the overall system is partitioned by means of two moving pistons, they may be constrained by the condition that $C_2^{\infty ,\mathrm{str}}\left(x\right)=V_1+V_2+V_3=V$ is constant or has an externally imposed time dependence. Another example is the set of stoichiometric relations between species concentrations which, in a detailed chemical or electrochemical kinetic model of a closed system, is imposed by the assumed scheme of reaction mechanisms.

In addition, for most systems of practical interest, such as for those of interest in thermodynamics and chemical kinetics, the system admits one or more functionals that we denote by $C_i^{\infty ,\mathrm{dyn}}\left(x\right)$, which are time-invariant under the internal dynamics $D^{\mathrm{int}}\left(x\right)$, i.e., such that for every state x,

$${\left.\frac{d{C}_i^{\infty ,\mathrm{dyn}}\left(x\right)}{dt}\right|}_{\mathrm{int}}=\sum _{j=1}^{n_{\mathrm{x}}}{a}_{ij}^{\infty ,\mathrm{dyn}}\left(x\right)D_j^{\mathrm{int}}\left(x\right)=0$$

(4)

where we defined the gradient vectors

$$a_{ij}^{\infty ,\mathrm{dyn}}\left(x\right)=\frac{\partial C_i^{\infty ,\mathrm{dyn}}\left(x\right)}{\partial x_j}$$

(5)

The functionals $C_i^{\infty ,\mathrm{dyn}}\left(x\right)$ are said to represent the constants of the motion of the system. For example, the energy and the number of atomic nuclei of each type are constants of the motion for a chemically reacting open system. They are time invariant under the internal dynamics, but their values may change as a result of external mechanisms such as those of work, heat, bulk flow, and diffusion interactions.

If there are r such constants of the motion, with linearly independent gradient vectors $a_{ij}^{\infty ,\mathrm{dyn}}\left(x\right)$, the r Equation (4) imply that the set of Equation (2) which define the equilibrium states are not all linearly independent and, therefore, the solution depends on at least r arbitrary parameters. In other words, the equilibrium states form at least an r-parameter family.

For each given state x, it is useful to identify, within the combined set of elemental constraints $C_i^{\infty ,\mathrm{str}}\left(x\right)$ and constants of the motion $C_i^{\infty ,\mathrm{dyn}}\left(x\right)$, a complete subset that we denote simply by $C_i^{\infty }\left(x\right)$ with $i=1,\cdots ,n_{\mathrm{c}}^{\infty }$, such that the $n_{\mathrm{c}}^{\infty }$ gradient vectors $a_{ij}^{\infty }\left(x\right)$ are linearly independent. The subset must be complete in the sense that any other functional in the full set has a gradient vector which is a linear combination of the vectors $a_{ij}^{\infty }\left(x\right)$.

For many systems of practical interest, there is a state functional $F\left(x\right)$ that is non-decreasing under the internal dynamics of the system [70], i.e., is such that

$${\left.\frac{dF\left(x\right)}{dt}\right|}_{\mathrm{int}}=\sum _{j=1}^{n_{\mathrm{x}}}{\Gamma }_j\left(x\right)D_j^{\mathrm{int}}\left(x\right)\ge 0$$

(6)

for every x, where we defined the gradient vector Γ with components

$${\Gamma }_j\left(x\right)=\frac{\partial F\left(x\right)}{\partial x_j}$$

(7)

Our next step is to assume that the functional $F\left(x\right)$ is such that in the neighborhood of each one of the equilibrium states of interest (the stable equilibrium states, in thermodynamics), the functional

$$\Delta F\left(x\right)=F\left(x^{\mathrm{eq}}\right)-F\left(x\right)$$

(8)

when restricted to the subset of states

$$\left\{x \mathrm{with} C_i^{\infty }\left(x\right)=C_i^{\infty }\left(x^{\mathrm{eq}}\right)\mathrm{for} i=1,\cdots ,n_{\mathrm{c}}^{\infty }\right\}$$

(9)

is a Lyapunov functional [71] and, hence, each equilibrium state is conditionally stable [71] with respect to perturbations that do not alter the values of the $C_i^{\infty }$’s.

Clearly, $\Delta F\left(x\right)$ can be a Lyapunov functional when restricted on the set (9) only if $\Delta F\left(x\right)\ge \alpha \left(\left|x-x^{\mathrm{eq}}\right|\right)$ where α is some strictly-increasing positive function with $\alpha \left(0\right)=0$ and, hence, the equilibrium state $x^{\mathrm{eq}}$ is the unique state which maximizes the functional $F\left(x\right)$ on the set Equation (9).

To fix ideas, by analogy with the thermodynamics context, we will assume that the functional $F\left(x\right)$ is twice differentiable and at least locally strictly concave in its variables, so that the Hessian matrix of generalized stiffness moduli

$$\frac{\partial {\Gamma }_j\left(x\right)}{\partial x_m}=\frac{{\partial }^{2}F\left(x\right)}{\partial x_m\partial x_j}$$

(10)

is well defined, negative semi-definite, and clearly symmetric by Schwarz’s theorem on mixed second derivatives. Moreover it has a nonzero determinant and, since it is also the Jacobian matrix of the system $\Gamma =\Gamma \left(x\right)$ of $n_{\mathrm{x}}$ relations

$${\Gamma }_j={\Gamma }_j\left(x\right)$$

(11)

it follows that the relations $\Gamma =\Gamma \left(x\right)$ can be inverted to yield the system $x=\widehat{x}(\Gamma )$ of $n_{\mathrm{x}}$ relations

$$x_m={\widehat{x}}_m(\Gamma )$$

(12)

Using this, we can express the functional $F\left(x\right)$ also as a functional of the vector Γ, i.e., $\widehat{F}(\Gamma )=F\left(\widehat{x}(\Gamma )\right)$, for which x is the gradient vector, with components

$${\widehat{x}}_m(\Gamma )=\frac{\partial \widehat{F}(\Gamma )}{\partial {\Gamma }_m}$$

(13)

and the symmetric, negative semi-definite Hessian matrix of generalized capacities,

$$\frac{\partial {\widehat{x}}_m(\Gamma )}{\partial {\Gamma }_j}=\frac{{\partial }^2\widehat{F}(\Gamma )}{\partial {\Gamma }_j\partial {\Gamma }_m}$$

(14)

is the inverse of the generalized stiffness matrix in Equation (10).

As a result of the above assumptions, the equilibrium states can be found by solving the following constrained maximization problem

$$\begin{array}{c}{x}^{\mathrm{eq}} \mathrm{maximizes} F\left(x\right) \mathrm{subject} \mathrm{to}\\ C_i^{\infty }\left(x\right)=c_i^{\infty }i=1,\dots ,n_{\mathrm{c}}^{\infty }\end{array}$$

(15)

where $c_i^{\infty }$ is either a fixed value of a geometrical constraint or the value of a constant of the motion. We have implicitly assumed that the functionals $C_i^{\infty }\left(x\right)$ form a complete set of geometrical constraints and constants of the motion with linearly independent gradients $a_{1j}^{\infty }\left(x\right),\dots ,a_{n_{\mathrm{c}}^{\infty }j}^{\infty }\left(x\right)$. The solution $x^{\mathrm{eq}}$ of the maximization problem satisfies the set of Euler-Lagrange equations

$${\Gamma }_j\left(x^{\mathrm{eq}}\right)=\sum _{i=1}^{n_{\mathrm{c}}^{\infty }}{\gamma }_i^{\infty }{a}_{ij}^{\infty }\left(x^{\mathrm{eq}}\right)$$

(16)

where ${\gamma }_i^{\infty }$ is the Lagrange multiplier associated with the i-th constraint. Together with the constraints (15), Equation (16) yields an $n_{\mathrm{c}}^{\infty }$-parameter family of equilibrium states

$$x^{\mathrm{eq}}=x^{\mathrm{eq}}(c_1^{\infty },\dots ,c_{n_{\mathrm{c}}^{\infty }}^{\infty })$$

(17)

which are stable with respect to all perturbations [71] if the function $x^{\mathrm{eq}}\left(c^{\infty }\right)$ is continuous in each $c_i^{\infty }$. As is well-known, the standard procedure to obtain (17) is as follows. First we invert Equation (16) to express the state variables in terms of the multipliers, $x^{\mathrm{eq}}=x^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)$; these relations are then inserted into (15) to express the values of the constraints as functions of the multipliers, $c^{\infty }=C^{\infty }\left(x^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)\right)$; these are then inverted to express the multipliers in terms of the values of the constraints, ${\mathbf{\gamma }}^{\infty }={\mathbf{\gamma }}^{\infty }\left(c^{\infty }\right)$, which inserted into $x^{\mathrm{eq}}=x^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)$ yields $x^{\mathrm{eq}}=x^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\left(c^{\infty }\right)\right)$, i.e., Equation (17). Moreover, it is noteworthy that as a result of the last inversion, we have the identity

$$\begin{array}{ccc}& {\displaystyle \sum _{m=1}^{n_{\mathrm{c}}^{\infty }}\frac{\partial C_i^{\infty }\left(x^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)\right)}{\partial {\gamma }_m^{\infty }}\frac{\partial {\gamma }_m^{\infty }\left(c^{\infty }\right)}{\partial c_k^{\infty }}=}& \hfill \end{array}$$

$$\begin{array}{ccc}& {\displaystyle \sum _{m=1}^{n_{\mathrm{c}}^{\infty }}\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{ij}^{\infty }\frac{\partial x_j^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)}{\partial {\gamma }_m^{\infty }}\frac{\partial {\gamma }_m^{\infty }\left(c^{\infty }\right)}{\partial c_k^{\infty }}={\delta }_{ik}}& \hfill \end{array}$$

(18)

which can also be viewed as a consequence of the condition that the constraints have linearly independent gradients, as necessary for the assumed invertibilities.

Using Equation (16) and Equation (18) we obtain the relations

$$\begin{array}{ccc}& {\displaystyle \frac{\partial F\left(x^{\mathrm{eq}}\left(c^{\infty }\right)\right)}{\partial c_k^{\infty }}=\sum _{j=1}^{n_{\mathrm{sp}}}{\Gamma }_j\frac{\partial x_j^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\left(c^{\infty }\right)\right)}{\partial c_k^{\infty }}=}& \hfill \end{array}$$

$$\begin{array}{ccc}& {\displaystyle \sum _{j=1}^{n_{\mathrm{sp}}}\sum _{i=1}^{n_{\mathrm{c}}^{\infty }}{\gamma }_i^{\infty }{a}_{ij}^{\infty }\sum _{m=1}^{n_{\mathrm{c}}^{\infty }}\frac{\partial x_j^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)}{\partial {\gamma }_m^{\infty }}\frac{\partial {\gamma }_m^{\infty }\left(c^{\infty }\right)}{\partial c_k^{\infty }}={\gamma }_k^{\infty }\left(c^{\infty }\right)}& \hfill \end{array}$$

(19)

which express, geometrically, the general orthogonality condition between the vector ${\mathbf{\gamma }}^{\infty }$ of Lagrange multipliers and the constant-F contour surfaces in $c^{\infty }$ space.

An immediate consequence of Equation (19) and Schwarz’s theorem on mixed second derivatives is the following Maxwell relations expressing the reciprocity of the cross dependence of the constraint-potentials on the values of the constraints,

$$\frac{\partial {\gamma }_i^{\infty }\left(c^{\infty }\right)}{\partial c_k^{\infty }}=\frac{{\partial }^{2}F\left(x^{\mathrm{eq}}\left(c^{\infty }\right)\right)}{\partial c_i^{\infty }\partial c_k^{\infty }}=\frac{\partial {\gamma }_k^{\infty }\left(c^{\infty }\right)}{\partial c_i^{\infty }}$$

(20)

Geometrically, the Hessian matrix of $F\left(x^{\mathrm{eq}}\left(c^{\infty }\right)\right)$ represents the curvature matrix of the F contour surfaces in $c^{\infty }$ space.

Similarly, the dual representation of the same surface in ${\mathbf{\gamma }}^{\infty }$ space is obtained by defining the Legendre transform of $F\left(x^{\mathrm{eq}}\left(c^{\infty }\right)\right)$ as follows,

$$G\left({\mathbf{\gamma }}^{\infty }\right)=\sum _{i=1}^{n_{\mathrm{c}}^{\infty }}{\gamma }_i^{\infty }{C}_i^{\infty }\left(x^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)\right)-F\left(x^{\mathrm{eq}}\left({\mathbf{\gamma }}^{\infty }\right)\right)$$

(21)

where the family of stable equilibrium states is now re-parametrized in terms of the $n_{\mathrm{c}}^{\infty }$ Lagrange multipliers,

$$x^{\mathrm{eq}}=x^{\mathrm{eq}}({\gamma }_1^{\infty },\dots ,{\gamma }_{n_{\mathrm{c}}^{\infty }}^{\infty })$$

(22)

It is easy to show that in ${\mathbf{\gamma }}^{\infty }$ space the orthogonality condition is between the vector $c^{\infty }$ of constraint values and the constant-G contour surfaces,

$$\frac{\partial G\left({\mathbf{\gamma }}^{\infty }\right)}{\partial {\gamma }_k^{\infty }}=c_k^{\infty }\left({\mathbf{\gamma }}^{\infty }\right)$$

(23)

Moreover, we have the dual Maxwell relations expressing the reciprocity of the cross dependence of the values of the constraints on the constraint-potentials,

$$\frac{\partial c_i^{\infty }\left(c^{\infty }\right)}{\partial {\gamma }_k^{\infty }}=\frac{{\partial }^{2}G\left({\mathbf{\gamma }}^{\infty }\right)}{\partial {\gamma }_i^{\infty }\partial {\gamma }_k^{\infty }}=\frac{\partial c_k^{\infty }\left(c^{\infty }\right)}{\partial {\gamma }_i^{\infty }}$$

(24)

and we can readily show that the Hessian matrix of G in Equation (24), which represents the curvature matrix of the G contour surfaces in ${\mathbf{\gamma }}^{\infty }$ space, is the inverse of the Hessian matrix of F in Equation (20).

If the internal dynamics is much faster than the external, we may assume that for all practical purposes the system always passes through equilibrium states, namely,

$$x\left(t\right)=x^{\mathrm{eq}}\left(c^{\infty }\left(t\right)\right)$$

(25)

where the values of the constants of the motion vary according to Equation (1), i.e., in view of Equation (4),

$${\dot{c}}_i^{\infty }\left(t\right)={\dot{C}}_i^{\infty }\left(x\left(t\right)\right)=\sum _{j=1}^{n_{\mathrm{x}}}{a}_{ij}^{\infty }\left(x\left(t\right)\right)D_j^{\mathrm{ext}}(x\left(t\right),t)$$

(26)

For example, this is a reasonable assumption for the local description of non-reacting flow fields subject to moderate spatial gradients. Actually, recent molecular dynamics simulations [72] have shown this to be an excellent approximation even for very large spatial gradients. But it would not be a reasonable assumption for reacting flows with reaction time scales of the same order as the transport time scales.

3. Constrained Equilibrium Method

In general, the evolution of the state, $x\left(t\right)$, is obtained by solving the differential Equation (1) for a given initial state $x\left(0\right)$ and the exogenous dynamics $D^{\mathrm{ext}}(x,t)$. However, the detailed modeling of $D^{\mathrm{int}}\left(x\right)$ may be a formidable task.

Upon gaining a deeper understanding of the dynamics of a system, it is often discovered that many of the governing kinetic mechanisms are of secondary importance, such that the time dependent behavior of the system is relatively insensitive to them. For example the three body reactions dominate the exothermic dynamics of the process of expansion of combustion products within a supersonic nozzle or in an internal combustion engine, while bimolecular reactions are almost equilibrated [35,73]. Also, the almost fuel independent burning velocity of the majority of the hydrocarbon fuels suggests the existence of an intermediate constrained-equilibrium path in which the fuel molecule breaks down into smaller fragments, mostly C1, C2 and C3, which could drastically simplify kinetic modeling of laminar flame propagation.

Ideally, by scrutinizing the features of the full endogenous dynamics of a system, we could classify the various governing mechanisms in a hierarchical structure based on their time scale as follows. We would single out all the internal mechanisms that act on a time scale shorter than a given scale τ, and identify a complete set of functionals, additional to the $C_i^{\infty }\left(x\right)$’s, that are also time invariant if the internal dynamics of the system is restricted to only those mechanisms that have a time scale shorter than τ. We denote these additional time-invariant functionals by

$$C_k^{\tau }\left(x\right)k=1,\dots ,n_{\mathrm{c}}^{\tau }$$

(27)

and we impose the condition that their gradient vectors,

$$a_{kj}^{\tau }\left(x\right)=\frac{\partial C_k^{\tau }\left(x\right)}{\partial x_j}$$

(28)

are linearly independent of the gradients $a_{kj}^{\infty }\left(x\right)$. Clearly, as τ increases, the number $n_{\mathrm{c}}^{\tau }$ of such functionals decreases to reach $n_{\mathrm{c}}^{\infty }=0$. Conversely, as τ decreases and reaches the limit $\tau =0$, at which time scale no mechanism can be active, the number $n_{\mathrm{c}}^{\tau }$ increases to reach $n_0=n_{\mathrm{x}}-n_{\mathrm{c}}^{\infty }$ where $n_{\mathrm{x}}$ is the number of state variables in vector x and $n_{\mathrm{c}}^{\infty }$ is the number of (τ-independent) time invariants $C_k^{\infty }\left(x\right)$.

In principle, one would have to examine a hierarchy of dynamical equations of the form

$$\dot{x}=D^{\mathrm{int},\le \tau }\left(x\right)$$

(29)

where $D^{\mathrm{int},\le \tau }\left(x\right)$ accounts for all the internal mechanisms with time scale shorter than τ. The functionals $C_k^{\tau }\left(x\right)$ are the additional constants of the motion for $D^{\mathrm{int},\le \tau }\left(x\right)$ with gradients $a_{kj}^{\tau }\left(x\right)$ that are linearly independent of the gradients $a_{ij}^{\infty }\left(x\right)$ of the constants of the motion of the full internal dynamics $D^{\mathrm{int},\le \infty }\left(x\right)=D^{\mathrm{int}}\left(x\right)$. This means that no mechanism with time scale shorter than τ is capable of altering the value of the functionals $C_k^{\tau }\left(x\right)$.

If we are interested in describing the dynamics of the system only on a time scale longer than a certain scale τ, then we are not interested in the details of the time evolution for time scales shorter than τ. We may therefore assume that the fast mechanisms, described by $D^{\mathrm{int},\le \tau }\left(x\right)$, contribute to the overall evolution by letting the functional $F\left(x\right)$ increase in practically no time to the highest value compatible with the instantaneous values of the functionals $C_i^{\infty }\left(x\right)$—which are constants of the motion for the overall dynamics—and the functionals $C_k^{\tau }\left(x\right)$—which are constants of the motion for the fast dynamics. Thus, we only need to develop a detailed model of the dynamical mechanisms acting on time scales between τ and ∞, which we will denote by the symbol $D^{>\tau }\left(x\right)$. In other words, we assume that for all practical purposes the system is at all times in some constrained-equilibrium state $x^{\mathrm{ce}}$ that can be found by solving the following constrained maximization problem

$$x^{\mathrm{ce}} \mathrm{maximizes} F\left(x\right) \mathrm{subject} \mathrm{to}$$

(30)

$$C_i^{\infty }\left(x\right)=c_i^{\infty }i=1,\dots ,n_{\mathrm{c}}^{\infty }$$

(31)

$$C_k^{\tau }\left(x\right)=c_k^{\tau }k=1,\dots ,n_{\mathrm{c}}^{\tau }$$

(32)

The solution $x^{\mathrm{ce}}$ of the maximization problem satisfies the set of Euler-Lagrange equations

$${\Gamma }_j\left(x^{\mathrm{ce}}\right)=\sum _{i=1}^{n_{\mathrm{c}}^{\infty }}{\gamma }_i^{\infty }{a}_{ij}^{\infty }\left(x^{\mathrm{ce}}\right)+\sum _{k=1}^{n_{\mathrm{c}}^{\tau }}{\gamma }_k^{\tau }{a}_{kj}^{\tau }\left(x^{\mathrm{ce}}\right))$$

(33)

which, together with the constraints (31) and (32), yields an $(n_{\mathrm{c}}^{\infty }+n_{\mathrm{c}}^{\tau })$-parameter family of constrained-equilibrium states

$$x^{\mathrm{ce}}=x^{\mathrm{ce}}(c_1^{\infty },\dots ,c_{n_{\mathrm{c}}^{\infty }}^{\infty },c_1^{\tau },\dots ,c_{n_{\mathrm{c}}^{\tau }}^{\tau })$$

(34)

Again, Equation (33) imply that, within this family, the function $F\left(x^{\mathrm{ce}}(c^{\infty },c^{\tau })\right)$ satisfies the following orthogonality relations

$$\frac{\partial F\left(x^{\mathrm{ce}}(c^{\infty },c^{\tau })\right)}{\partial c_i^{\infty }}={\gamma }_i^{\infty },\frac{\partial F\left(x^{\mathrm{ce}}(c^{\infty },c^{\tau })\right)}{\partial c_k^{\tau }}={\gamma }_k^{\tau }$$

(35)

where ${\gamma }_i^{\infty }$ and ${\gamma }_k^{\tau }$ are the constraint-potentials conjugated to the constraint functionals $C_i^{\infty }\left(x\right)$ and $C_k^{\tau }\left(x\right)$. Therefore, the relations and geometrical representations discussed in the previous section extend here as well for the constant-F surfaces of the constrained equilibrium states in $(c^{\infty },c^{\tau })$ space and for the dual constant-G surfaces in $({\mathbf{\gamma }}^{\infty },{\mathbf{\gamma }}^{\tau })$ space. As a result of this duality, the time dependence can be followed both in the space of the constraints and in that of the constraint-potentials.

3.1. Rate Equations for the Constraints

If the exogenous dynamics, $D^{\mathrm{ext}}(x,t)$, operates on a time scale larger than τ, then we see that the evolution of the system can be described by modeling only those internal mechanisms with time scale longer than τ. Indeed,

$$\begin{array}{ccc}& {\displaystyle {\dot{c}}_i^{\infty }\left(t\right)=\sum _{j=1}^{n_{\mathrm{x}}}{a}_{ij}^{\infty }\left(x\left(t\right)\right)D_j^{\mathrm{ext}}(x\left(t\right),t)}& \end{array}$$

(36)

$$\begin{array}{cc}& {\displaystyle {\dot{c}}_k^{\tau }\left(t\right)=\sum _{j=1}^{n_{\mathrm{x}}}{a}_{kj}^{\tau }\left(x\left(t\right)\right)\left[D_j^{\mathrm{int},>\tau }\left(x\left(t\right)\right)+D_j^{\mathrm{ext}}(x\left(t\right),t)\right]}\end{array}$$

(37)

To follow the evolution of the state variable in $(c^{\infty },c^{\tau })$ space, we assume at each instant in time

$$x\left(t\right)=x^{\mathrm{ce}}(c_1^{\infty }\left(t\right),\dots ,c_{n_{\mathrm{c}}^{\infty }}^{\infty }\left(t\right),c_1^{\tau }\left(t\right),\dots ,c_{n_{\mathrm{c}}^{\tau }}^{\tau }\left(t\right))$$

(38)

where the function $x^{\mathrm{ce}}(c^{\infty },c^{\tau })$ is that found in Equation (34) as a result of the maximization just discussed.

As a result, the solution of this scheme is based on the set of $n_{\mathrm{c}}^{\infty }+n_{\mathrm{c}}^{\tau }$ explicit differential Equations (36) and (37) coupled with the $n_{\mathrm{c}}^{\infty }+n_{\mathrm{c}}^{\tau }+n_{\mathrm{x}}$ Equations (31)–(33), where the $n_{\mathrm{x}}+2n_{\mathrm{c}}^{\infty }+2n_{\mathrm{c}}^{\tau }$ unknown functions of time are $x_j\left(t\right)$, $c_i^{\infty }\left(t\right)$, $c_k^{\tau }\left(t\right)$, ${\gamma }_i^{\infty }\left(t\right)$, and ${\gamma }_k^{\tau }\left(t\right)$.

In particular, Equation (37) shows that the only internal mechanisms that determine the rates of change of the values of the functionals $C_i^{\tau }\left(x\right)$ are those acting on time scales between τ and ∞. Thus, the set of Equations (36)–(38) is independent of the details of the dynamical mechanisms whose time scale is shorter than the given τ.

It is finally noteworthy that Relations (4) and (6) impose important restrictions on the expressions of $D^{\mathrm{int},>\tau }\left(x\right)$ that can be adopted to model the rate-controlling slow endogenous dynamics of the system. The first condition that we have already used in Equation (36) is that for every constrained-equilibrium state $x^{\mathrm{ce}}$,

$${\dot{c}}_i^{\infty }{|}_{\mathrm{int}}=\sum _{j=1}^{n_{\mathrm{x}}}{a}_{ij}^{\infty }\left(x^{\mathrm{ce}}\right)D_j^{\mathrm{int},>\tau }\left(x^{\mathrm{ce}}\right)=0i=1,\dots ,n_{\mathrm{c}}^{\infty }$$

(39)

namely, the functionals $C_i^{\infty }\left(x\right)$ must be time invariants of the endogenous dynamics also in the constrained-equilibrium scheme. The second, more demanding condition is that the state functional $F\left(x\right)$ must be non-decreasing under the endogenous dynamics. Setting the exogenous term $D^{\mathrm{ext}}(x\left(t\right),t)=0$ in Equation (36) and (37) and using Relations (35), this condition requires that at any constrained-equilibrium state $x^{\mathrm{ce}}\left(t\right)$ along the assumed time evolution,

$$\begin{array}{ccc}\hfill {\left.\frac{dF\left(x^{\mathrm{ce}}\right)}{dt}\right|}_{\mathrm{int}}& =& \sum _{j=1}^{n_{\mathrm{x}}}{\Gamma }_j\left(x^{\mathrm{ce}}\right)D_j^{\mathrm{int},>\tau }\left(x^{\mathrm{ce}}\right)\hfill \end{array}$$

$$\begin{array}{ccc}& =& \sum _{k=1}^{n_{\mathrm{c}}^{\tau }}{\gamma }_k^{\tau }\left(x^{\mathrm{ce}}\right)\sum _{j=1}^{n_{\mathrm{x}}}{a}_{kj}^{\tau }\left(x^{\mathrm{ce}}\right)D_j^{\mathrm{int},>\tau }\left(x^{\mathrm{ce}}\right)\hfill \end{array}$$

$$\begin{array}{ccc}& =& \sum _{k=1}^{n_{\mathrm{c}}^{\tau }}{\gamma }_k^{\tau }\left(x^{\mathrm{ce}}\right){\dot{c}}_k^{\tau }{|}_{\mathrm{int}}\left(x^{\mathrm{ce}}\right)\ge 0\hfill \end{array}$$

(40)

where in writing the second equality we used Equation (33) and Equation (39), and in the third we used Equation (37).

For example, like in the application of RCCE to chemical kinetics that we discuss in the next section, a fully consistent assumption is to use, in Equation (37), $D_j^{\mathrm{int},>\tau }\left(x^{\mathrm{ce}}\right)=D_j^{\mathrm{int}}\left(x^{\mathrm{ce}}\right)$, i.e., to use the full internal dynamics (the detailed kinetic scheme) to evaluate the rates of change of the slow constraints once the constrained equilibrium state $x^{\mathrm{ce}}$ is known. Then, the nonnegativity of Equation (40) is guaranteed in general by the assumed nonnegativity of Equation (6). See also note [62].

3.2. Rate Equations for the Constraint-Potentials

Although direct integration of the rate-equations (36) and (37) for the constraints is relatively straightforward and simple to implement, it has proved to be relatively inefficient and time consuming if, at each time step during the integration, the constrained-equilibrium values (36) of the state variables corresponding to the values of the constraints are recalculated by solving Equation (33) using a general purpose maximization code that does not accept as initial guess the composition computed at the previous time step [52]. For the purpose of an implementation of the method, also an alternative method has been proposed [53] and implemented [34,50,54,74] which, instead of the direct integration of the rate-equations for the constraint-potentials, solves the dual problem in terms of an implicit set of rate equations for the constraint-potentials ${\gamma }_i$. A general advantage of this method is the reduction by $n_{\mathrm{c}}$ (from $n_{\mathrm{c}}+n_{\mathrm{x}}$ to $n_{\mathrm{x}}$) of the number of implicit equations that need to be solved during the integration.

For simplicity, let us rewrite the constraints without separating them into two sets, so that (31) and (32) are united in the single set of constraints,

$$C_i\left(x^{\mathrm{ce}}\right)=c_{i}i=1,\dots ,n_{\mathrm{c}}$$

(41)

and the Euler-Lagrange equations are

$${\Gamma }_j\left(x^{\mathrm{ce}}\right)=\sum _{k=1}^{n_{\mathrm{c}}}{\gamma }_k{a}_{kj}\left(x^{\mathrm{ce}}\right)$$

(42)

and implicitly define the relation

$$x^{\mathrm{ce}}=x^{\mathrm{ce}}\left(\mathbf{\gamma }\right)$$

(43)

By differentiating Equation (42) with respect to time and using (43) in the right-hand side, we may write

$$\sum _{n=1}^{n_{\mathrm{x}}}\frac{\partial {\Gamma }_m}{\partial x_n}{\dot{x}}_n^{\mathrm{ce}}=\sum _{k=1}^{n_{\mathrm{c}}}{\dot{\gamma }}_k{a}_{km}+\sum _{i=1}^{n_{\mathrm{c}}}{\gamma }_i\sum _{n=1}^{n_{\mathrm{x}}}\frac{\partial a_{im}}{\partial x_n}\sum _{k=1}^{n_{\mathrm{c}}}\frac{\partial x_n^{\mathrm{ce}}}{\partial {\gamma }_k}{\dot{\gamma }}_k$$

(44)

Defining for convenience the matrix

$$a_{km}^{\prime }=a_{km}+\sum _{i=1}^{n_{\mathrm{c}}}{\gamma }_i\sum _{n=1}^{n_{\mathrm{x}}}\frac{\partial a_{im}}{\partial x_n}\frac{\partial x_n^{\mathrm{ce}}}{\partial {\gamma }_k}$$

(45)

multiplying both sides of (44) by $\partial {\widehat{x}}_j/\partial {\Gamma }_m$, summing over m and recalling that relations (11) and (12) are one the inverse of the other, so that ${\sum }_m(\partial {\widehat{x}}_j/\partial {\Gamma }_m)(\partial {\Gamma }_m/\partial x_n)={\delta }_{jn}$, we obtain

$${\dot{x}}_j^{\mathrm{ce}}=\sum _{m=1}^{n_{\mathrm{x}}}\frac{\partial {\widehat{x}}_j}{\partial {\Gamma }_m}\sum _{k=1}^{n_{\mathrm{c}}}{\dot{\gamma }}_k{a}_{km}^{\prime }$$

(46)

We now use Equation (1) with x given by Equation (43) to evaluate ${\dot{x}}_j^{\mathrm{ce}}$ in Equation (46), we multiply both sides by $a_{ij}$, sum over j, and after a minor rearrangement and re-inserting all dependences explicitly, we obtain the rate equations

$$\sum _{k=1}^{n_{\mathrm{c}}}{\dot{\gamma }}_k\sum _{j=1}^{n_{\mathrm{x}}}\sum _{m=1}^{n_{\mathrm{x}}}{a}_{ij}\left(x^{\mathrm{ce}}\left(\mathbf{\gamma }\right)\right){\left.\frac{\partial {\widehat{x}}_j(\Gamma )}{\partial {\Gamma }_m}\right|}_{\Gamma =\mathbf{\gamma }·a\left(x^{\mathrm{ce}}\left(\mathbf{\gamma }\right)\right)}{a}_{km}^{\prime }\left(x^{\mathrm{ce}}\left(\mathbf{\gamma }\right)\right)$$

$$=\sum _{j=1}^{n_{\mathrm{x}}}{a}_{ij}\left(x^{\mathrm{ce}}\left(\mathbf{\gamma }\right)\right)\left[D_j^{\mathrm{int}}\left(x^{\mathrm{ce}}\left(\mathbf{\gamma }\right)\right)+D_j^{\mathrm{ext}}(x^{\mathrm{ce}}\left(\mathbf{\gamma }\right),t)\right]$$

(47)

which form a set of $n_{\mathrm{c}}$ differential equations for the $n_{\mathrm{c}}$ constraint potentials $\mathbf{\gamma }$.

In summary, the comparison between the dual solution methods presented in this and the preceding section may be clarified by writing the equations to be solved in compact vector notation. The method based on the rate equations for the constraints requires the solution of a set of $n_{\mathrm{c}}$ explicit differential equations coupled with $n_{\mathrm{c}}$ explicit equations plus $n_{\mathrm{x}}$ implicit (usually transcendental) equations, with the forms

$$\begin{array}{ccc}\hfill \dot{c}\left(t\right)& =& a\left(x\right(t\left)\right)·D\left(x\right(t\left)\right)\hfill \end{array}$$

(48)

$$\begin{array}{ccc}\hfill c\left(t\right)& =& C\left(x\right(t\left)\right)\hfill \end{array}$$

(49)

$$\begin{array}{ccc}\hfill \Gamma \left(x\right(t\left)\right)& =& \mathbf{\gamma }\left(t\right)·a\left(x\right(t\left)\right)\hfill \end{array}$$

(50)

for the unknown functions $c\left(t\right)$, $\mathbf{\gamma }\left(t\right)$, and $x\left(t\right)$.

The method based on the rate equations for the constraint-potentials requires the solution of a set of $n_{\mathrm{c}}$ implicit differential equations coupled with $n_{\mathrm{x}}$ implicit (usually transcendental) equations, with the forms

$$\begin{array}{ccc}\hfill d\left(x\left(t\right)\right)·\dot{\mathbf{\gamma }}\left(t\right)& =& a\left(x\right(t\left)\right)·D\left(x\right(t\left)\right)\hfill \end{array}$$

(51)

$$\begin{array}{ccc}\hfill \Gamma \left(x\right(t\left)\right)& =& \mathbf{\gamma }\left(t\right)·a\left(x\right(t\left)\right)\hfill \end{array}$$

(52)

for the unknown functions $\mathbf{\gamma }\left(t\right)$ and $x\left(t\right)$.

4. Rate-Controlled Constrained-Equilibrium in Chemical Kinetics of Gas Mixtures

The development of models for describing the dynamic evolution of chemically reacting systems is a fundamental objective of chemical kinetics. As explained in the introduction, the most accurate way of doing this involves: (i) specifying the state and species variables to be included in the model; (ii) developing a detailed kinetic model (DKM) that is comprehensive, hierarchical, and includes all the important reaction pathways governing the state variables; (iii) obtaining accurate rate data for each individual reaction in the DKM either using experimental techniques or using state-of-the-art quantum chemistry methods [75]; (iv) compiling a “full set” of rate-equations for these variables based on the postulated DKM; and (v) integrating this set of equations to obtain the time-dependent behavior of the system.

For hydrocarbon oxidation, the problem is that the DKMs involving C/H/O/N molecules can easily involve thousands of chemical species and isomers, and tens of thousands of chemical reactions for heavy molecules. Aside from the expensive computational effort required to treat such models in realistic systems involving reacting turbulent flows, one major modeling challenge is to assign reliable rate data to every individual reaction, especially the intermediate reactions involving heavy radicals and branched molecules, for which even the current most accurate quantum chemistry calculations offer factor-of-ten accuracies, and virtually many of which are irrelevant after all.

It is anticipated that the RCCE method, which was originally developed and is mainly understood and applied as a model order reduction technique, could significantly reduce the modeling effort required in developing a DKM. One hidden, not very well received promising potential of RCCE is the ability of bottom-up modeling. One starts with an exhaustive list of species and elemental constraints, which are sufficient to fix the equilibrium state. No kinetic mechanism is required at this stage in the RCCE model development, whereas the minimum number of linearly independent reactions in the corresponding DKM is equal to the number of species only to fix the equilibrium state. The appreciable gain at this very first stage of modeling is more pronounced when heavy molecules including several thousands of species are considered, formation of PAH and soot for example. More constraints and more reactions can be added to improve the accuracy to the desired level. If the only constraints are those imposed by slowly changing state variables, the RCCE method is equivalent to a “shifting” Local Thermodynamic (complete chemical) Equilibrium (LTE) calculation. If the number of constraints in an RCCE model is identical to the number of species in a DKM model, then they are, for all practical purposes, identical.

RCCE fundamentals will be summarized below by starting from the standard treatment of chemical reactions, briefly reviewed in Section 4.1 to Section 4.6, and then by focusing on the RCCE formulation for a homogeneous system (Section 4.7) and on two examples (Section 5).

4.1. Equilibrium Properties of Non-Reacting Gas Mixtures

In homogeneous, closed-reactor chemical kinetics or RCCE, where the usual frozen chemical equilibrium assumption is made to calculate the thermodynamic properties, the Lyapunov functional is the stable-equilibrium fundamental relation for the entropy S of a non-reacting system with the instantaneous values of the internal energy U, the volume V, and the numbers $N_j$ of particles of each different species,

$$F\left(x\right)=\tilde{S}\left(x\right)=\frac{1}{R}S(U,V,N_1,\cdots ,N_{n_{\mathrm{sp}}})$$

(53)

where $R=N_{\mathrm{Av}}{k}_{\mathrm{B}}$ with $N_{\mathrm{Av}}$ Avogadro’s number and $k_{\mathrm{B}}$ Boltzmann’s constant. Therefore, the state vector is $x=(U,V,N_1,\cdots ,N_{n_{\mathrm{sp}}})$. In this context, the components of the gradient vector Γ are the potentials conjugated with the state variables,

$$\Gamma =\left[\frac{1}{RT},\frac{p}{RT},-\frac{{\mu }_1}{RT},\cdots ,-\frac{{\mu }_{n_{\mathrm{sp}}}}{RT}\right]$$

(54)

where T is the temperature, p the pressure, and ${\mu }_j$ the chemical potential or partial molar Gibbs free energy of species j in the mixture. Under the typical simple-system approximation ([36], p. 263) whereby the function $S(U,V,N_1,\cdots ,N_{n_{\mathrm{sp}}})$ is homogeneous of first degree in all its variables, the validity of the Euler relation, $U=TS-pV+{\sum }_j{\mu }_j{N}_j$, entails a relation among the potentials, which are, therefore, not all independent of one another. For this reason, it is often more convenient to select specific properties as state variables, i.e., depending on the context, either

• on a molar basis with

  $$F^N\left(x\right)=\frac{\tilde{S}\left(x\right)}{N}=\tilde{s}(u,v,x_1,\cdots ,x_{n_{\mathrm{sp}}-1})$$

  (55)

  and $\tilde{s}=\tilde{S}/N$, $u=U/N$, $v=V/N$, $x_j=N_j/N$, $N={\sum }_j{N}_j$, so that in this context the state vector is $x=(u,v,x_1,\cdots ,x_{n_{\mathrm{sp}}-1})$ and it is easy to show that the gradient vector is

  $${\Gamma }^N=\left[\frac{1}{RT},-\frac{p}{RT},-\frac{{\mu }_1-{\mu }_{n_{\mathrm{sp}}}}{RT},\cdots ,-\frac{{\mu }_{n_{\mathrm{sp}}-1}-{\mu }_{n_{\mathrm{sp}}}}{RT}\right]$$

  (56)
• on a mass basis with

  $$F^m\left(x\right)=\frac{\tilde{S}\left(x\right)}{m}={\tilde{s}}^m(u^m,\rho ,y_1,\cdots ,y_{n_{\mathrm{sp}}-1})$$

  (57)

  and ${\tilde{s}}^m=\tilde{S}/m$, $u^m=U/m$, $\rho =m/V$ is the density, $y_j=M_j{N}_j/m$ and $M_j$ the mass fraction and molar mass of species j, and $m={\sum }_j{N}_j{M}_j$ the mass, so that in this context the state vector is $x=(u^m,\rho ,y_1,\cdots ,y_{n_{\mathrm{sp}}-1})$ and it is easy to show that the gradient vector is

  $${\Gamma }^m=\left[\frac{1}{RT},-\frac{p}{{\rho }^{2}RT},-\frac{{\mu }_1^m-{\mu }_{n_{\mathrm{sp}}}^m}{RT},\cdots ,-\frac{{\mu }_{n_{\mathrm{sp}}-1}^m-{\mu }_{n_{\mathrm{sp}}}^m}{RT}\right]$$

  (58)

  where ${\mu }_j^m={\mu }_j/M_j$; or
• on a volume basis with

  $$F^V\left(x\right)=\frac{\tilde{S}\left(x\right)}{V}={\tilde{s}}^V(u^V,\left[N_1\right],\cdots ,\left[N_{n_{\mathrm{sp}}}\right])$$

  (59)

  and ${\tilde{s}}^V=S/V$, $u^V=U/V$, $\left[N_j\right]=N_j/V$ is the concentration of species j, so that in this context the state vector is $x=(u^V,\left[N_1\right],\cdots ,\left[N_{n_{\mathrm{sp}}}\right])$, and it is easy to show that the gradient vector is

  $${\Gamma }^V=\left[\frac{1}{RT},-\frac{{\mu }_1}{RT},\cdots ,-\frac{{\mu }_{n_{\mathrm{sp}}}}{RT}\right]$$

  (60)

  This latter case is convenient in non-homogeneous, reacting and non-reacting flows, where fluid mechanics and chemical kinetics or RCCE are combined with the methods of non-equilibrium thermodynamics. The pressure is obtained from the Euler relation, $p/RT={\tilde{s}}^V-u^V/RT+{\sum }_{j=1}^{n_{\mathrm{sp}}}\left[N_j\right]{\mu }_j/RT$.

We start, therefore, with a system consisting of $n_{\mathrm{sp}}$ different species confined in a volume V, with numbers of particles, denoted by $N_j$, sufficiently large so that the simple system approximation applies.

For a non-reacting mixture in the absence of external fields, the stable-equilibrium-state fundamental relation for a simple system with volume as the only parameter has the form

$$S=S_{\mathrm{off}}(U,V,N_1,\cdots ,N_{n_{\mathrm{sp}}})$$

(61)

where the subscript “off” is a reminder [76] that the functional relation is that of a non-reacting mixture (all reactions “turned off”).

The Gibbs relation, which is the differential of (61),

$$dS=\frac{1}{T}dU+\frac{p}{T}dV-\sum _{j=1}^{n_{\mathrm{sp}}}\frac{{\mu }_j}{T}d{N}_j$$

(62)

allows us to write the rate equation which relates the rates of change of entropy, energy, volume, and composition when the system is modeled as evolving along a continuous time sequence of states, each of which would be a stable equilibrium state if the all reactions were inhibited. Introducing, for convenience, the “dimensionless entropy” $\tilde{S}$, the “entropic temperature” β, and the dimensionless “entropic chemical potential” ${\lambda }_j$ of species j,

$$\tilde{S}=\frac{S}{R},\beta =\frac{1}{RT},{\lambda }_j=-\frac{{\mu }_j}{RT}$$

(63)

and assuming that U, V, and the $N_j$’s are functions of time, and at any instant in time the entropy is given by Equation (61), it is easy to see that Equation (62) implies the following relation among the time rates of change,

$$\frac{d\tilde{S}}{dt}=\beta \frac{dU}{dt}+\beta p\frac{dV}{dt}+\sum _{j=1}^{n_{\mathrm{sp}}}{\lambda }_j\frac{d{N}_j}{dt}$$

(64)

When applied to an infinitesimal control volume in a flow field, the set of assumptions which lead to Equation (64) constitutes the so-called “local equilibrium” assumption of fluid mechanics. In this and the following sections, instead, we restrict the treatment to the simpler case of a homogeneous system. But it is important to note that the “shifting equilibrium” assumption is adopted also when we deal with reacting mixtures, for which the thermodynamic properties are assumed to be those of a non-reacting mixture at stable (frozen) equilibrium with the same energy, volume, and composition.

4.2. Balance Equations for Species, Energy, and Entropy

Let us now write the rate equations expressing the species, energy, and entropy balances, using the following general notation

$$\begin{array}{ccc}\hfill \frac{d{N}_j}{dt}& =& -{\dot{N}}_j^{\to }+V{\dot{\omega }}_j\hfill \end{array}$$

(65)

$$\begin{array}{ccc}\hfill \frac{dU}{dt}& =& -{\dot{E}}^{\to }-p\frac{dV}{dt}+V\Phi \hfill \end{array}$$

(66)

$$\begin{array}{ccc}\hfill \frac{dS}{dt}& =& -{\dot{S}}^{\to }+V\sigma \hfill \end{array}$$

(67)

where ${\dot{N}}_j^{\to }$, ${\dot{E}}^{\to }$, and ${\dot{S}}^{\to }$ represent the net rates of outflow (inflow if negative) of species, energy, and entropy due to interactions across the system’s boundaries such as heat, bulk flow, diffusion, and work other than the expansion (or compression) work (which is accounted for explicitly in view of the “shifting equilibrium” assumption). Moreover, ${\dot{\omega }}_j$, Φ, and σ represent the net rates per unit volume of “generation within the system” of particles of type j, internal energy, and entropy due, respectively, to chemical reactions, viscous dissipation, and irreversibilities.

Substitution of Equations (65)–(67) in Equation (64) yields the following expression for the rate of entropy generation by irreversibility per unit volume,

$$\sigma ={\sigma }_{\mathrm{chem}}+{\sigma }_{\mathrm{flux}}+{\sigma }_{\mathrm{diss}}$$

(68)

where

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{chem}}& =& R\sum _{j=1}^{n_{\mathrm{sp}}}{\lambda }_j{\dot{\omega }}_j\ge 0\hfill \end{array}$$

(69)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{flux}}& =& R\left[{\dot{\tilde{S}}}^{\to }-\beta {\dot{E}}^{\to }-\sum _{j=1}^{n_{\mathrm{sp}}}{\lambda }_j{\dot{N}}_j^{\to }\right]\ge 0\hfill \end{array}$$

(70)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{diss}}& =& R\beta \Phi \ge 0\hfill \end{array}$$

(71)

and, of course, ${\dot{\tilde{S}}}^{\to }={\dot{S}}^{\to }/R$.

The principle of entropy non-decrease, which follows from the Second Law, imposes that${\sigma }_{\mathrm{chem}}+{\sigma }_{\mathrm{flux}}+{\sigma }_{\mathrm{diss}}\ge 0$. The reason why we write the stronger condition that ${\sigma }_{\mathrm{chem}}$, ${\sigma }_{\mathrm{flux}}$, and ${\sigma }_{\mathrm{diss}}$ must be independently non-negative, is to be found in the Curie principle of nonequilibrium thermodynamics, whereby chemical reactions, which are scalar processes, cannot couple to heat and mass fluxes, which are vectorial, nor to momentum transfer, which is tensorial.

4.3. Complete Chemical Equilibrium

If the composition is allowed to vary subject only to conservation of the atomic nuclei present in the initial composition $N_{j,\mathrm{initial}}$, thereby including the formation mechanisms of every possible chemical species that can be conceivably assembled with the given atomic nuclei, then the corresponding stable-equilibrium composition is the so-called “complete chemical equilibrium state”. For a given value of the energy U, the volume V, and the initial numbers of particles $N_{j,\mathrm{initial}}$ of the different species, the complete chemical equilibrium state is found by maximizing the entropy subject to the following constraints, which express conservation of the numbers of atomic nuclei of each kind that of course no chemical reaction mechanism can change,

$$C_i^{\infty }\left(x\right)=\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{ij}^{\infty }{N}_j=c_i^{\infty }fori=1,\cdots ,n_{\mathrm{el}}$$

(72)

where $a_{ij}^{\infty }$ is the number of nuclei of type i in one molecule of species j, $n_{\mathrm{el}}$ is the number of elemental species which can be formed with the given initial set of particles, and the superscript ∞ is used here because we will see below that in RCCE these parameters appear in elemental constraints of the model.

For a closed system, the value of each $c_i^{\infty }$ is determined once and for all by the “initial” composition $N_{j,\mathrm{initial}}$, and it cannot be altered by the chemistry.

In terms of the terminology of Section 2 and the examples discussed at the beginning of Section 4.1, to proceed without making further modeling assumptions, it is convenient to adopt the state variables $x=(U,V,N_1,\cdots ,N_{n_{\mathrm{sp}}})$. Therefore, the complete chemical equilibrium state is obtained by maximizing the function $F\left(x\right)=\tilde{S}\left(x\right)$ as given by Equation (53) subject to given values of U and V and the constraint Equation (72). Introducing multipliers ${\gamma }_U^{\infty }$, ${\gamma }_V^{\infty }$, and ${\gamma }_i^{\infty }$ for each constraint, and maximizing the Lagrangian

$$L\left(x\right)=\tilde{S}-{\gamma }_V^{\infty }V-{\gamma }_U^{\infty }U-\sum _{i=1}^{n_{\mathrm{el}}}{\gamma }_i^{\infty }{C}_i^{\infty }\left(x\right)$$

(73)

yields the necessary conditions

$$\frac{\partial L}{\partial U}=\beta -{\gamma }_U^{\infty }=0 \mathrm{that} \mathrm{is} {\gamma }_U^{\infty }=\beta$$

(74)

$$\frac{\partial L}{\partial V}=\beta p-{\gamma }_V^{\infty }=0 \mathrm{that} \mathrm{is} {\gamma }_V^{\infty }=\beta p$$

(75)

$$\frac{\partial L}{\partial N_j}={\lambda }_j-\sum _{i=1}^{n_{\mathrm{el}}}{\gamma }_i^{\infty }{a}_{ij}^{\infty }=0$$

(76)

Because in general ${\lambda }_j={\lambda }_j(U,V,N_1,\cdots ,N_{n_{\mathrm{sp}}})$, the latter conditions imply the relations

$${\lambda }_j(U,V,N_1,\cdots ,N_{n_{\mathrm{sp}}})=\sum _{i=1}^{n_{\mathrm{el}}}{\gamma }_i^{\infty }{a}_{ij}^{\infty }$$

(77)

which can be inverted at fixed U and V to give

$$N_j=N_j(U,V,{\gamma }_1,\cdots ,{\gamma }_{n_{\mathrm{el}}})$$

(78)

When these are inserted in the constraints (72), the resulting relations can be inverted at fixed U and V to yield

$${\gamma }_i={\gamma }_i(U,V,c_1^{\infty },\cdots ,c_{n_{\mathrm{el}}}^{\infty })$$

(79)

which, in turn, inserted into Equation (78) yields the complete chemical composition corresponding to the given values of U, V, and initial number of nuclei $c^{\infty }$,

$$N_j=N_j(U,V,c_1^{\infty },\cdots ,c_{n_{\mathrm{el}}}^{\infty })$$

(80)

So far, no specific assumption has been made regarding the relations between the stable-equilibrium properties of the mixture and those of the component species when pure. This is where further assumptions, for example to model a particular non-ideal mixture behavior valid in dense gas, liquid phase, or phase change, should be considered as additional “elemental” constraints of the model. In the next section, we exemplify this by assuming an ideal mixture of ideal gases, so as to recover the standard treatment of chemical kinetics and the original formulation of RCCE.

It is noteworthy that because of conditions (74) and (75), the term $\tilde{S}-\beta U-\beta pV$ in the right-hand side of the Lagrangian (73) coincides with $-G/RT$ where G is the Gibbs free energy. Therefore, maximizing $\tilde{S}$ subject to given values of U and V is equivalent to minimizing the dimensionless Planck function $\tilde{K}=G/RT$. Thus, again, for non-ideal mixture behavior, we must proceed by adopting one of the many model expressions available in the literature (see, e.g., [77,78]) for the Gibbs free energy G or the Helmholtz free energy $F=G-pV$ as functions of the other state variables.

4.4. Complete Chemical Equilibrium for Ideal Gas Mixtures

Let us now make and discuss the following standard assumptions in the treatment of ideal gas mixtures: (i) simple system model; (ii) Gibbs-Dalton ideal mixture behavior; (iii) ideal gas behavior for each of the component species when pure and at the same temperature of the mixture, in the entire range between the standard reference pressure $p_o$ and the pressure of the mixture.

The simple system model ([36], p. 266) implies the validity of the Euler relation which can be written either as $U=TS-pV+{\sum }_j{\mu }_j{N}_j$ or as

$$\tilde{S}=\beta U+\beta pV+\sum _{j=1}^{n_{\mathrm{sp}}}{\lambda }_j{N}_j$$

(81)

The Gibbs–Dalton ideal mixture model ([36], p. 481) implies the following relations for the mixture properties S, U, V, and the pure substance properties of the component species,

$$\begin{array}{ccc}\hfill S& =& \sum _{j=1}^{n_{\mathrm{sp}}}{N}_j{s}_{jj}(\beta ,p_j)\hfill \end{array}$$

(82)

$$\begin{array}{ccc}\hfill U& =& \sum _{j=1}^{n_{\mathrm{sp}}}{N}_j{u}_{jj}(\beta ,p_j)\hfill \end{array}$$

(83)

$$\begin{array}{ccc}\hfill V& =& N_j{v}_{jj}(\beta ,p_j)foreveryj=1,\cdots ,n_{\mathrm{sp}}\hfill \end{array}$$

(84)

where $p_j$ denotes the partial pressure of species j in the mixture and $s_{jj}$, $u_{jj}$, and $v_{jj}$ the mole specific properties of pure species j, where the double-j subscript is a reminder that these properties refer to the species j when pure. Under this model, the Dalton law applies, $p={\sum }_{j=1}^{n_{\mathrm{sp}}}{p}_j$.

Finally, assuming ideal gas behavior, $v_{jj}(\beta ,p)=\beta p$, for every species j and all values of β and p, then the partial pressures are given by $p_j=N_j/\beta V$, and Equations (82)–(84) may be written in the form

$$\begin{array}{ccc}& {\displaystyle \tilde{S}=\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left[{\tilde{s}}_{jj}(\beta ,p_o)-ln\frac{p}{p_o}-ln\frac{N_j}{N}\right]}& \end{array}$$

(85)

$$\begin{array}{ccc}& {\displaystyle C_U^{\infty }\left(x\right)=\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j{u}_{jj}\left(\beta \right)=U}& \end{array}$$

(86)

$$\begin{array}{cc}& {\displaystyle C_V^{\infty }\left(x\right)=\frac{1}{\beta p}\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j=V}\end{array}$$

(87)

where $p_o$ is the standard reference pressure. We further note that, in our notation, the relations ${\mu }_{jj}=h_{jj}-T{s}_{jj}$, $h_{jj}=u_{jj}+RT$ translate into

$${\tilde{s}}_{jj}={\lambda }_{jj}+\beta u_{jj}+1$$

(88)

so that we may rewrite Equation (85) as

$$\tilde{S}=\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left[{\lambda }_{jj}(\beta ,p_o)+\beta u_{jj}\left(\beta \right)+1-ln\frac{p}{p_o}-ln\frac{N_j}{N}\right]$$

(89)

As a result, a natural set of state variables in this context is

$$x=(\beta ,p,N_1,\cdots ,N_{n_{\mathrm{sp}}})$$

(90)

and so the complete chemical equilibrium state is obtained by maximizing the function $F\left(x\right)=\tilde{S}\left(x\right)$ as given by Equation (85) subject to the constraints given by Equations (72), (86), and (87). Introducing multipliers ${\gamma }_i^{\infty }$, ${\gamma }_U^{\infty }$, and ${\gamma }_V^{\infty }$ for these constraints, and maximizing the Lagrangian

$$L\left(x\right)=\tilde{S}\left(x\right)-{\gamma }_V^{\infty }{C}_V^{\infty }\left(x\right)-{\gamma }_U^{\infty }{C}_U^{\infty }\left(x\right)-\sum _{i=1}^{n_{\mathrm{el}}}{\gamma }_i^{\infty }{C}_i^{\infty }\left(x\right)$$

(91)

yields the necessary conditions

$$\frac{\partial L}{\partial p}=N\left[-\frac{1}{p}+{\gamma }_V^{\infty }\frac{1}{\beta p^2}\right]=0 \mathrm{that} \mathrm{is} {\gamma }_V^{\infty }=\beta p$$

(92)

$$\frac{\partial L}{\partial \beta }=\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left[-\frac{{\tilde{c}}_{pjj}}{\beta }+{\gamma }_U^{\infty }\frac{{\tilde{c}}_{vjj}}{{\beta }^2}+\frac{1}{\beta }\right]=0 \mathrm{that} \mathrm{is} {\gamma }_U^{\infty }=\beta$$

(93)

$$\frac{\partial L}{\partial N_j}={\lambda }_{jj}(\beta ,p_o)-ln\frac{p}{p_o}-ln\frac{N_j}{N}-\sum _{i=1}^{n_{\mathrm{el}}}{\gamma }_i^{\infty }{a}_{ij}^{\infty }=0$$

(94)

where we made use of the relations $\partial {\tilde{s}}_{jj}/\partial \beta =-{\tilde{c}}_{pjj}/\beta$, $\partial u_{jj}/\partial \beta =-{\tilde{c}}_{vjj}/{\beta }^2$, ${\tilde{c}}_{pjj}={\tilde{c}}_{vjj}+1$, where of course ${\tilde{c}}_p=c_p/R$ and ${\tilde{c}}_v=c_v/R$. Using Equation (87) to eliminate N in Equation (94), finally yields the relation

$$N_j=\beta p_{o}Vexp\left[{\lambda }_{jj}(\beta ,p_o)-\sum _{i=1}^{n_{\mathrm{el}}}{\gamma }_i^{\infty }{a}_{ij}^{\infty }\right]$$

(95)

Equations (72), (86), (87), and (95) form a system of $n_{\mathrm{sp}}+n_{\mathrm{el}}+2$ implicit equations which, for given values of U, V, and the $c_i^{\infty }$’s, can be solved to yield the complete chemical equilibrium values of the $n_{\mathrm{sp}}+2$ state variables β, p, $N_j$, and the $n_{\mathrm{el}}$ constraint-potentials ${\gamma }_i^{\infty }$.

If the values of the energy U, the volume V, and the number of atomic nuclei of each type, $c_i^{\infty }$, are allowed to vary in time, the system will evolve through a sequence of complete-equilibrium states at a rate controlled by the imposed rates. This is the well-known case of “shifting complete chemical equilibrium”, already denoted by LTE.

For completeness, it is convenient to recall and write in more standard notation, that each (dimensionless) ${\lambda }_{jj}(\beta ,p_o)$ in Equation (95) is computed from the molar Gibbs free energies of formation of the j-th species at standard conditions and the temperature dependence of its molar specific heat capacity at constant pressure $c_{pjj}\left(T\right)$, via the relation

$${\lambda }_{jj}(1/RT,p_o)=-\frac{\Delta g_{\mathrm{f}}^o{)}_{jj}}{RT}-\frac{1}{RT}\left[{\int }_{T_o}^T{c}_{pjj}\left(T\right)dT-T{\int }_{T_o}^T\frac{c_{pjj}\left(T\right)}{T}dT\right]$$

(96)

Similarly each molar specific internal energy $u_{jj}\left(\beta \right)$ may be computed via the relation

$$u_{jj}(1/RT)={(\Delta h_{\mathrm{f}}^o)}_{jj}-RT+{\int }_{T_o}^T{c}_{pjj}\left(T\right)dT$$

(97)

4.5. Chemical Equilibrium of Ideal Gas Mixtures with Respect to a Set of Reactions

If the composition is allowed to vary according to one or more chemical reaction mechanisms, the maximization is to be done for specified values of U and V and the allowed stoichiometric variations in composition. So, let us assume that the DKM is based on the following $n_r$ reaction mechanisms,

$$\sum _{j=1}^{n_{\mathrm{sp}}}{\nu }_{j\ell }^{+}{B}_j\leftrightarrow \sum _{j=1}^{n_{\mathrm{sp}}}{\nu }_{j\ell }^{-}{B}_{j}for\ell =1,\cdots ,n_r$$

(98)

where $B_j$ denotes the chemical symbol of the j-th species, and the net stoichiometric coefficients ${\nu }_{j\ell }={\nu }_{j\ell }^{-}-{\nu }_{j\ell }^{+}$ satisfy the atomic nuclei conservation identities

$$\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{ij}^{\infty }{\nu }_{j\ell }=0\mathrm{for} \mathrm{each} \mathrm{element} i \mathrm{and} \mathrm{reaction} \ell$$

(99)

where again $a_{ij}^{\infty }$ is the number of nuclei of type i in one molecule of species j.

As a result of this assumption, the term $V{\dot{\omega }}_j$ in the species balance equation (65), which expresses the rate of formation (consumption, if negative) of particles of species j due to the allowed chemical reaction mechanisms, must satisfy the so-called proportionality relations

$$V{\dot{\omega }}_j=\sum _{\ell =1}^{n_{\mathrm{r}}}{\nu }_{j\ell }{\dot{\epsilon }}_{\ell }$$

(100)

where ${\epsilon }_{\ell }$ is the so-called coordinate or degree of advancement of the ℓ-th reaction.

The Second Law requires that at stable equilibrium the entropy be maximal within each set of states that have the same values of the energy U and the volume V, and compositions that are compatible (through the allowed reaction mechanisms) with the same initial composition $N_{j,\mathrm{initial}}$. The set of compositions compatible with $N_{j,\mathrm{initial}}$ is the r-parameter family of compositions found by integrating in time Equation (65) with ${\dot{N}}_j^{\to }=0$ and with $V{\dot{\omega }}_j$ given by Equation (100). This integration, if we set for convenience ${\epsilon }_{\ell }\left(0\right)$, yields

$$N_j\left(\mathbf{\epsilon }\right)=N_{j,\mathrm{initial}}+\sum _{\ell =1}^{n_{\mathrm{r}}}{\nu }_{j\ell }{\epsilon }_{\ell }$$

(101)

Equation (101) together with the additional trivial restriction $N_j\ge 0$ defines through the r parameters ${\epsilon }_{\ell }$ the family of compatible compositions over which the entropy maximization is to be restricted. Inserting this restriction into Equations (85)–(87), they become,

$$\begin{array}{ccc}& {\displaystyle \tilde{S}=\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left(\mathbf{\epsilon }\right)\left[{\tilde{s}}_{jj}(\beta ,p_o)-ln\frac{p}{p_o}-ln\frac{N_j\left(\mathbf{\epsilon }\right)}{N\left(\mathbf{\epsilon }\right)}\right]}& \end{array}$$

(102)

$$\begin{array}{ccc}& {\displaystyle C_U^{\infty }\left(x\right)=\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left(\mathbf{\epsilon }\right)u_{jj}\left(\beta \right)=U}& \end{array}$$

(103)

$$\begin{array}{cc}& {\displaystyle C_V^{\infty }\left(x\right)=\frac{1}{\beta p}\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left(\mathbf{\epsilon }\right)=V}\end{array}$$

(104)

Therefore, a natural set of state variables in this more constrained context is

$$x=(\beta ,p,{\epsilon }_1,\cdots ,{\epsilon }_{n_{\mathrm{r}}})$$

(105)

and so the chemical equilibrium state is obtained by maximizing the function $F\left(x\right)=\tilde{S}\left(x\right)$ as given by Equation (102) subject only to the constraints (103) and (104). In fact, we note that the element-conservation constraint (72) is automatically satisfied because the stoichiometric coefficients satisfy Equation (99).

Introducing multipliers ${\gamma }_U^{\infty }$ and ${\gamma }_V^{\infty }$ for the two constraints, and maximizing the Lagrangian

$$L\left(x\right)=\tilde{S}\left(x\right)-{\gamma }_V^{\infty }{C}_V^{\infty }\left(x\right)-{\gamma }_U^{\infty }{C}_U^{\infty }\left(x\right)$$

(106)

yields the necessary conditions ${\gamma }_V^{\infty }=\beta p$, ${\gamma }_U^{\infty }=\beta$ as in Equation (92) and Equation (93), and

$$\frac{\partial L}{\partial {\epsilon }_{\ell }}=\sum _{j=1}^{n_{\mathrm{sp}}}{\nu }_{j\ell }\left[{\lambda }_{jj}(\beta ,p_o)-ln\frac{p}{p_o}-ln\frac{N_j}{N}\right]=0$$

(107)

where we used $\partial N_j/\partial {\epsilon }_{\ell }={\nu }_{j\ell }$ which follows from Equation (101).

By defining the dimensionless equilibrium “constant” of reaction ℓ at entropic temperature β,

$$K_{\ell }^o\left(\beta \right)=exp\left[\sum _{j=1}^{n_{\mathrm{sp}}}{\lambda }_{jj}(\beta ,p_o){\nu }_{j\ell }\right]$$

(108)

and ${\nu }_{\ell }={\sum }_{j=1}^{n_{\mathrm{sp}}}{\nu }_{j\ell }$, Equation (107) takes the familiar form in terms of the species concentrations $\left[N_j\right]=N_j/V$, namely, at equilibrium,

$$\prod _{j=1}^{n_{\mathrm{sp}}}{\left[N_j\right]}^{{\nu }_{j\ell }}={\left(\beta p_o\right)}^{{\nu }_{\ell }}{K}_{\ell }^o\left(\beta \right)$$

(109)

which is the well-known law of mass action. The right-hand side of Equation (109) defines the so-called equilibrium constant based on concentrations,

$$K_{\ell }^{\mathrm{c}o}\left(\beta \right)={\left(\beta p_o\right)}^{{\nu }_{\ell }}{K}_{\ell }^o\left(\beta \right)$$

(110)

whose dimensions depend on the value of ${\nu }_{\ell }$.

We finally note that the potential ${\Gamma }_{\ell }$ conjugated to the reaction coordinate ${\epsilon }_{\ell }$,

$${\Gamma }_{\ell }=\frac{\partial \tilde{S}\left(x\right)}{\partial {\epsilon }_{\ell }}=\sum _{j=1}^{n_{\mathrm{sp}}}{\nu }_{j\ell }{\lambda }_j$$

(111)

is the dimensionless form of the entropic affinity of reaction ℓ, i.e., with the notation of [76], ${\Gamma }_{\ell }=Y_{\ell }/R$ and with the usual notation ${\Gamma }_{\ell }={\mathcal{A}}_{\ell }/RT$ where ${\mathcal{A}}_{\ell }=-{\sum }_{j=1}^{n_{\mathrm{sp}}}{\nu }_{j\ell }{\mu }_j$ is the de Donder affinity of the reaction.

For the discussion in the next section, it is important to note that ${\Gamma }_{\ell }$ may be expressed, using Equation (102), as

$${\Gamma }_{\ell }=ln{K}_{\ell }^o\left(\beta \right)-ln\frac{{\prod }_{j=1}^{n_{\mathrm{sp}}}{\left[N_j\right]}^{{\nu }_{j\ell }}}{{\left(\beta p_o\right)}^{{\nu }_{\ell }}}$$

(112)

Therefore, the chemical equilibrium condition (109) is equivalently expressed by the condition ${\Gamma }_{\ell }=0$, i.e., that the affinity of every reaction vanishes at equilibrium.

However, the entropic affinity ${\Gamma }_{\ell }$ is well defined also off chemical equilibrium. In fact, inserting Equations (100) in Equation (69) and using Equation (111), the rate of entropy generation due to the r allowed chemical reactions, may be written as

$${\sigma }_{\mathrm{chem}}=R\sum _{j=1}^{n_{\mathrm{sp}}}{\lambda }_j{\dot{\omega }}_j=\frac{R}{V}\sum _{\ell =1}^{n_{\mathrm{r}}}{\Gamma }_{\ell }{\dot{\epsilon }}_{\ell }$$

(113)

In the next section, we show that under the standard kinetic modeling assumptions, the Second Law requirement that ${\sigma }_{\mathrm{chem}}$ be positive semi-definite implies that ${\Gamma }_{\ell }$ can be interpreted as a (dimensionless) measure of the degree of disequilibrium with respect to the ℓ-th reaction.

4.6. Entropy Production and Detailed Balance

The rate of change of reaction coordinate ${\epsilon }_{\ell }$ per unit volume is usually called the net rate of the ℓ-th reaction, denoted by $r_{\ell }$ and expressed as the difference between the forward and reverse reaction rates,

$$r_{\ell }={\dot{\epsilon }}_{\ell }/V=r_{\ell }^{+}-r_{\ell }^{-}$$

(114)

So, the entropy generation rate density is also given by

$${\sigma }_{\mathrm{chem}}=R\sum _{\ell =1}^{n_{\mathrm{r}}}(r_{\ell }^{+}-r_{\ell }^{-}){\Gamma }_{\ell }$$

(115)

For ideal gas mixture behavior, the standard model in chemical kinetics is based on the following expressions for the forward and reverse (positive) reaction rates

$$\begin{array}{ccc}\hfill r_{\ell }^{+}& =& k_{\ell }^{+}\prod _{j=1}^{n_{\mathrm{sp}}}{\left[N_j\right]}^{{\nu }_{j\ell }^{+}}\hfill \end{array}$$

(116)

$$\begin{array}{ccc}\hfill r_{\ell }^{-}& =& k_{\ell }^{-}\prod _{j=1}^{n_{\mathrm{sp}}}{\left[N_j\right]}^{{\nu }_{j\ell }^{-}}\hfill \end{array}$$

(117)

where $k_{\ell }^{+}$ and $k_{\ell }^{-}$ are the so-called (positive) rate constants of the reaction.

Recalling that ${\nu }_{j\ell }={\nu }_{j\ell }^{-}-{\nu }_{j\ell }^{+}$ and taking the ratio of Equations (117) and (116) yields the relation

$$\prod _{j=1}^{n_{\mathrm{sp}}}{\left[N_j\right]}^{{\nu }_{j\ell }}=\frac{r_{\ell }^{-}}{r_{\ell }^{+}}\frac{k_{\ell }^{+}}{k_{\ell }^{-}}$$

(118)

which is valid both at and off equilibrium. Inserting this relation in Equation (112), we obtain, in general,

$${\Gamma }_{\ell }=ln\frac{r_{\ell }^{+}}{r_{\ell }^{-}}+ln\left(K_{\ell }^{\mathrm{c}o}\left(\beta \right)\frac{k_{\ell }^{-}}{k_{\ell }^{+}}\right)$$

(119)

which inserted in Equation (115) yields, in general,

$${\sigma }_{\mathrm{chem}}=R\sum _{\ell =1}^{n_{\mathrm{r}}}(r_{\ell }^{+}-r_{\ell }^{-})\left[ln\frac{r_{\ell }^{+}}{r_{\ell }^{-}}+ln\left(K_{\ell }^{\mathrm{c}o}\left(\beta \right)\frac{k_{\ell }^{-}}{k_{\ell }^{+}}\right)\right]$$

(120)

Let us rewrite the expressions just obtained, in terms of two useful parameters. The first, suggested by one of us [53] as a measure of the departure of reaction ℓ from equilibrium, is the “degree of disequilibrium”

$${\varphi }_{\ell }=ln\frac{r_{\ell }^{+}}{r_{\ell }^{-}}$$

(121)

The second parameter, which we will see can be interpreted as the “degree of departure from detailed balance”, is

$${\Delta }_{\ell }=ln\left(K_{\ell }^{\mathrm{c}o}\left(\beta \right)\frac{k_{\ell }^{-}}{k_{\ell }^{+}}\right)$$

(122)

Clearly, Equation (119) and Equation (120) rewrite as

$${\Gamma }_{\ell }={\varphi }_{\ell }+{\Delta }_{\ell }$$

(123)

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{chem}}& =& R\sum _{\ell =1}^{n_{\mathrm{r}}}{r}_{\ell }^{+}[1-exp(-{\varphi }_{\ell })]({\varphi }_{\ell }+{\Delta }_{\ell })\hfill \end{array}$$

(124)

$$\begin{array}{ccc}& =& R\sum _{\ell =1}^{n_{\mathrm{r}}}{r}_{\ell }^{+}[1-exp(-{\Gamma }_{\ell }+{\Delta }_{\ell })]{\Gamma }_{\ell }\hfill \end{array}$$

(125)

So far, we made no assumption about the functional dependence of the rate constants $k_{\ell }^{+}$ and $k_{\ell }^{-}$, and hence of ${\Delta }_{\ell }$, on the state variables.

Close to chemical equilibrium, in the limit of vanishing net reaction rates $r_{\ell }=r_{\ell }^{+}-r_{\ell }^{-}$, i.e., if ${\varphi }_{\ell }\to 0$ for every ℓ, the affinities must approach their vanishing equilibrium values, i.e., ${\Gamma }_{\ell }\to 0$ for every ℓ, therefore Equation (123) implies ${\Delta }_{\ell }\to 0$, which is the so-called rate quotient law or detailed balance condition, relating the ratio of the forward and reverse rate constants for every reaction to the corresponding equilibrium constant,

$$\underset{{\varphi }_{\ell }\to 0}{lim}\frac{k_{\ell }^{+}}{k_{\ell }^{-}}=K_{\ell }^{\mathrm{c}o}\left(\beta \right)$$

(126)

In general, away from chemical equilibrium, by Equation (122) and Equation (123) we have

$$\frac{k_{\ell }^{+}}{k_{\ell }^{-}}=K_{\ell }^{\mathrm{c}o}\left(\beta \right)exp({\varphi }_{\ell }-{\Gamma }_{\ell })$$

(127)

and the question of the validity of the detailed balance condition ${\Delta }_{\ell }=0$ can be addressed in various ways [79], mainly depending on the modeling assumptions about the functional dependence of the rate constants $k_{\ell }^{+}$ and $k_{\ell }^{-}$. However, the most general and restrictive condition for the thermodynamic consistency of any set of modeling assumptions is that the expression (124) for ${\sigma }_{\mathrm{chem}}$ be positive semi-definite over the entire range of compatible compositions, no matter how far from chemical equilibrium.

For example, let us assume that the rate constants exhibit neither cross effects nor any dependence on the affinities, so that $\partial {\Delta }_{\ell }/\partial {\Gamma }_m=0$ for every ℓ and m. Then it follows that ${\Delta }_{\ell }=0$ is a necessary and sufficient condition for ${\sigma }_{\mathrm{chem}}$ to be positive semi-definite. That the condition is sufficient is a straightforward consequence of the nonnegativity of the function $[1-exp(-\Gamma \left)\right]\Gamma$. One way to prove that the condition is also necessary is to require that the minimum value of ${\sigma }_{\mathrm{chem}}$ given by (125) with respect to variations of Γ be nonnegative (it actually must be zero, since at equilibrium the entropy production must be zero). Denoting by $\overline{\Gamma }$ the location of the minimum, where $\partial {\sigma }_{\mathrm{chem}}/\partial {\Gamma }_{\ell }=0$ for every ℓ, it is easy to show that

$${\sigma }_{\mathrm{chem},\min }=-R\sum _{\ell =1}^{n_{\mathrm{r}}}{r}_{\ell }^{+}{\overline{\Gamma }}_{\ell }^{2}exp(-{\overline{\Gamma }}_{\ell }+{\overline{\Delta }}_{\ell })$$

(128)

For this to be nonnegative, we must have ${\overline{\Gamma }}_{\ell }=0$ for every ℓ, which inserted back into the condition for the minimum, $\partial {\sigma }_{\mathrm{chem}}/\partial {\Gamma }_{\ell }=0$, yields $1-exp\left({\Delta }_{\ell }\right)=0$ for every ℓ, that is, ${\Delta }_{\ell }=0$. An equivalent more direct proof obtains by differentiating Equation (123) with respect to ${\Gamma }_{\ell }$ under the assumed independence of ${\Delta }_{\ell }$ on the affinities. This yields $\partial {\varphi }_{\ell }/\partial {\Gamma }_{\ell }=1$ which, together with the condition that at chemical equilibrium both ${\Gamma }_{\ell }=0$ and ${\varphi }_{\ell }=0$, implies ${\Gamma }_{\ell }={\varphi }_{\ell }$ and, hence, ${\Delta }_{\ell }=0$.

As a result we find that the detailed balance condition,

$$\frac{k_{\ell }^{+}}{k_{\ell }^{-}}=K_{\ell }^{\mathrm{c}o}\left(\beta \right)$$

(129)

holds also off equilibrium. This means that if the forward rate constant $k_{\ell }^{+}$ is assumed, as typical, to be a function of temperature and pressure only, then Equation (129) implies that also the reverse rate constant $k_{\ell }^{-}$ must be a function of temperature and pressure only, with $k_{\ell }^{+}/k_{\ell }^{-}$ independent of pressure.

Conversely, if we assume that the forward and reverse rate constants $k_{\ell }^{+}$ and $k_{\ell }^{-}$ are functions of temperature and pressure only, with $k_{\ell }^{+}/k_{\ell }^{-}$ independent of pressure, i.e., we assume that no matter how far from the chemical equilibrium composition, $k_{\ell }^{+}/k_{\ell }^{-}$ is a function of β only, then because of Condition (126) such function must be equal to $K_{\ell }^{\mathrm{c}o}\left(\beta \right)$ implying that ${\Delta }_{\ell }=0$ at all compositions.

In summary, any set of assumptions about the dependences of the rate constants that implies the general validity of detailed balance off-equilibrium implies also that the expressions for the rate of entropy generation per unit volume reduce to the following equivalent positive semi-definite forms

$$\begin{array}{ccc}\hfill {\sigma }_{\mathrm{chem}}& =& R\sum _{\ell =1}^{n_{\mathrm{r}}}(r_{\ell }^{+}-r_{\ell }^{-})ln\frac{r_{\ell }^{+}}{r_{\ell }^{-}}\hfill \end{array}$$

(130)

$$\begin{array}{ccc}& =& R\sum _{\ell =1}^{n_{\mathrm{r}}}{r}_{\ell }^{+}[1-exp(-{\varphi }_{\ell })]{\varphi }_{\ell }\hfill \end{array}$$

(131)

$$\begin{array}{ccc}& =& R\sum _{\ell =1}^{n_{\mathrm{r}}}{r}_{\ell }^{+}[1-exp(-{\Gamma }_{\ell })]{\Gamma }_{\ell }\hfill \end{array}$$

(132)

We finally note that sometimes the term “detailed balance” is used with a different meaning, to express the condition $r_{\ell }^{-}=r_{\ell }^{+}$ of equality of the forward and reverse reaction rates, i.e., in our language, to the absence of disequilibrium, ${\varphi }_{\ell }=0$. However, it is clear that close enough to equilibrium $r_{\ell }^{-}\approx r_{\ell }^{+}$.

4.7. Rate-Controlled Constrained-Equilibrium Formulation for a Closed System

The Rate-Controlled Constrained Equilibrium (RCCE) model assumes a time evolution through a sequence of states each maximizing the entropy subject not only to the constraints that define the complete chemical equilibrium state (Section 4.3 and Section 4.4), but also to a set of additional kinetically-controlled constraints. The number of additional constraints may be large or small.

At one extreme is the case that we call Detailed Model (DM) with the largest possible number of additional constraints. These are the number of particles of a selection of $n_{\mathrm{sp}}-n_{\mathrm{el}}$ species among those considered in the underlying DKM scheme. The additional time-dependent constraints are $C_1^{\tau }\left(x\right)=N_1$,..., $C_{n_{\mathrm{sp}}-n_{\mathrm{el}}}^{\tau }\left(x\right)=N_{n_{\mathrm{sp}}-n_{\mathrm{el}}}$, so that the evolution of the system is controlled by the rate-equations for the species, plus of course the prescribed time changes of U and V. DM and DKM differ in that once we have the DM calculation, then we may at once calculate the amounts of an augmented set of species, including any other species which can be made from the given initial nuclei in the system, even if not included in the underlying DKM scheme.

However, the RCCE method is most interesting when it effectively reproduces DKM results with a number of time-dependent constraints that is of course higher than in LTE but much smaller than in DM. So far, RCCE calculations have been based on identifying slow varying functionals $C_k^{\tau }\left(x\right)$ that are linear in the species number of moles $N_j$, and therefore can be written in the form

$$C_k^{\tau }\left(x\right)=\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{kj}^{\tau }{N}_j=c_k^{\tau }\mathrm{for} k=1,\cdots ,n_{\mathrm{c}}^{\tau }$$

(133)

where $a_{kj}^{\tau }$ represents the marginal impact of a change in $N_j$ onto changing the value of the slowly varying functional. The rate of change of the value $c_k^{\tau }$ is obtained, for the closed system, by using the rates of change of the amounts $N_j$ which result from using, in the rate-equations (100), the reaction rates $r_{\ell }$ known from the DKM scheme.

Because they have identical linear structure, we may formally combine Equation (72) and Equation (133), into a single set of $n_{\mathrm{c}}=n_{\mathrm{el}}+n_{\mathrm{c}}^{\tau }$ constraint equations,

$$C_i\left(x\right)=\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{ij}{N}_j=c_i\mathrm{for} i=1,\cdots ,n_{\mathrm{c}}$$

(134)

So, the maximization procedure that determines the concentrations at time t is formally identical to that outlined in Section 4.3 except that the number of constraints here is higher, $n_{\mathrm{c}}$ instead of $n_{\mathrm{el}}$. The resulting concentrations are given therefore by the expressions

$$N_j=\beta p_{o}Vexp\left[{\lambda }_{jj}(\beta ,p_o)-\sum _{i=1}^{n_{\mathrm{c}}}{\gamma }_i{a}_{ij}\right]$$

(135)

which differ from Equation (95) only in the absence of the superscript ∞, in accordance with the new, larger set of constraints (134).

Equations (89), (87), (134), and (135) form a system of $n_{\mathrm{sp}}+n_{\mathrm{c}}+2$ implicit equations which can be solved for given values of U, V, and the $c_i$’s, to yield the constrained equilibrium values of the $n_{\mathrm{sp}}+2$ state variables β, p, $N_j$, and the $n_{\mathrm{c}}=n_{\mathrm{el}}+n_{\mathrm{c}}^{\tau }$ constraint-potentials ${\gamma }_i$.

As the values of the energy U, the volume V, the number of atomic nuclei of each type, $c_i^{\infty }$, and the slowly varying constraints $c_i^{\tau }$ are allowed to vary in time, the system will evolve through a sequence of constrained-equilibrium states at a rate controlled by the imposed rates.

In view of their particular functional dependences and because the constraint gradients $a_{ij}$ are constants, the RCCE equations can be integrated more efficiently by rewriting them as rate equations for the constraint-potentials and making use of the species and energy balance Equations (65) and (66), and the kinetic Equations (100), (114), (116), and (117) that determine the rate of change of the values $c_i$ of the constraints. By differentiating Equations (86), (87), (134), and (135) with respect to time, after few rearrangements corresponding to the general procedure outlined in Section 3.2 for the general case, we obtain the following set of rate equations

$$\sum _{k=1}^{n_{\mathrm{c}}}{\dot{\gamma }}_k\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{kj}{N}_j{a}_{ij}=-\frac{\dot{\beta }}{\beta }\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\beta u_{jj}{a}_{ij}+\frac{\dot{V}}{V}\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j{a}_{ij}+\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{ij}{\dot{N}}_j^{\to }-\left\{\begin{array}{c}0,\mathrm{for} i=1,\cdots ,n_{\mathrm{el}},\hfill \\ \\ {\dot{c}}_{i,\mathrm{chem}}^{\tau }(\beta ,p,N),\hfill \\ \mathrm{for} i=n_{\mathrm{el}}+1,\cdots ,n_{\mathrm{c}},\hfill \end{array}\right.$$

(136)

where

$${\dot{c}}_{i,\mathrm{chem}}^{\tau }(\beta ,p,N)=V\sum _{\ell =1}^{n_{\mathrm{r}}}\left(\sum _{m=1}^{n_{\mathrm{sp}}}{a}_{im}{\nu }_{m\ell }\right)\left[\frac{k_{\ell }^{+}(\beta ,p)}{V^{{\nu }_{\ell }^{+}}}\prod _{j=1}^{n_{\mathrm{sp}}}{\left(N_j\right)}^{{\nu }_{j\ell }^{+}}-\frac{k_{\ell }^{-}(\beta ,p)}{V^{{\nu }_{\ell }^{-}}}\prod _{j=1}^{n_{\mathrm{sp}}}{\left(N_j\right)}^{{\nu }_{j\ell }^{-}}\right]$$

(137)

$$\sum _{k=1}^{n_{\mathrm{c}}}{\dot{\gamma }}_k\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{kj}{N}_j\beta u_{jj}=-\frac{\dot{\beta }}{\beta }\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left({\beta }^2{u}_{jj}^2+\frac{c_{vjj}}{R}\right)+\frac{\dot{V}}{V}\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\beta u_{jj}+\beta {\dot{E}}^{\to }+\beta p\dot{V}-\beta V\Phi$$

(138)

$$\sum _{k=1}^{n_{\mathrm{c}}}{\dot{\gamma }}_k\sum _{j=1}^{n_{\mathrm{sp}}}{a}_{kj}{N}_j=-\frac{\dot{\beta }}{\beta }\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j\left(\beta u_{jj}+1\right)-\frac{\dot{p}}{p}\sum _{j=1}^{n_{\mathrm{sp}}}{N}_j$$

(139)

These $n_{\mathrm{nc}}+2$ implicit differential equations together with the $n_{\mathrm{x}}$ Equation (135) can be solved for given values of ${\dot{E}}^{\to }$, $V\left(t\right)$, and the ${\dot{N}}_j^{\to }$, to yield the $n_{\mathrm{x}}+2$ state variables $\beta \left(t\right)$, $p\left(t\right)$, and $N_j\left(t\right)$, and the $n_{\mathrm{c}}$ constraint-potentials ${\gamma }_i\left(t\right)$.

It is finally important to note that the thermodynamic-consistency condition (40), that the entropy generation rate be non-negative definite under the endogenous dynamics, is always satisfied, because in evaluating the rates ${\dot{c}}_i^{\tau }$ we use the full internal dynamics (the detailed kinetic scheme), with forward and reverse rate constants $k_{\ell }^{+}$ and $k_{\ell }^{-}$ that depend on temperature and pressure only, a condition that we have seen in Section 4.6 warrants the nonnegativity of the entropy generation rate. As previously noted (see also [53,76]), whether or not there are kinetically-controlled constraints in addition to element conservation, all conceivable reactions contribute to the rate of entropy production. However, in the sum in Equation (120), the contribution of the non-constraint-changing reactions is through the adjustments in the chemical affinities of the constraint-changing reactions.

5. Selection of Constraints in Combustion Examples

The careful selection of constraints is the key to the success of the RCCE method. In general, the constraints must (i) include the elemental constraints; (ii) be linearly independent combinations of the species molar concentrations; (iii) constrain the system in the specified initial state; (iv) constrain global reactions in which reactants or intermediates go directly to products; and (v) determine the energy and entropy of the system within experimental accuracy. In addition, they should reflect whatever information is available about rate-limiting reactions which control the evolution of the system on the time scale of interest.

Under broad conditions in the gas phase, two important structural constraints are the total number of moles and the amount of free valence in a reacting system. The former is controlled by slow three-body dissociation/recombination reactions and the latter by reactions which make and break valence bonds. In connection with the ideas presented in Section 2 and Section 3, we observe that the total number of moles is a constant of the endogenous dynamics within the class of bimolecular reactions. Three-body reactions are generally slow in the endothermic direction because of the high activation energies required, and in the exothermic direction because of small three-body collision rates and small radical concentrations. They impose slowly varying constraints on the total molar concentration N and the free valence of the system. We denote these two constraints by the symbol M and FV, respectively. A finite value of FV is a necessary condition for chemical reactions to proceed.

A third important structural constraint, accounting for slow O−O bond-breaking reactions, is the free oxygen, FO, defined as any oxygen atom not directly bound to another oxygen atom. An increase in FO is a necessary condition for the formation of the major reaction products of hydrocarbon oxidation, H${}_2$O, CO${}_2$ and CO. Two additional constraints which slightly improve the agreement between RCCE and DKM calculations under some conditions are: OHO ≡ OH + O and DCO ≡ HCO + CO. The OHO constraint is a consequence of the relatively slow constraint-changing reaction RH + OH ↔ H${}_2$O + R coupled with the fast reaction RH + O ↔ OH + R which equilibrates OH and O. The DCO constraint is a consequence of the slow spin-forbidden reaction CO + HO${}_2$↔ CO${}_2$ + OH coupled with the fast reaction HCO + O${}_2$↔ CO + HO${}_2$ which equilibrates HCO and CO.

For systems involving the elements C, H and O, the five nontrivial constraints M, FV, FO, OHO, and DCO are independent of the initial reactants and may, therefore, be considered structurally “universal” constraints. Along with the three equilibrium reactions, H${}_2$ + O ↔ OH + H, H${}_2$ + HOO ↔ H${}_2$O${}_2$ + H, and HCO + O${}_2$↔ CO + HO${}_2$, they are sufficient to determine the constrained-equilibrium mole fractions of the 11 major hydrocarbon combustion products (H, O, OH, HO${}_2$, H${}_2$, O${}_2$, H${}_2$O, H${}_2$O${}_2$, HCO, CO and CO${}_2$) under both high and low temperature conditions.

For hydrocarbon oxidation, four additional fuel-dependent constraints have proved to be relevant [34]. The first is a constraint on the fuel, FU, imposed by slow hydrogen-abstraction reactions of the type FU + O${}_2$↔ FR + HO${}_2$ and even slower dissociation/recombination of the type AB + M ↔ A + B + M. This constraint is necessary to hold the system in its initial state. The second is a constraint on fuel radicals, FR, which is necessary to prevent the equilibration of forbidden global reactions such as CH${}_3$ + 2O${}_2$ + 2H${}_2$O ↔ CO${}_2$ + 2H${}_2$O${}_2$ + H${}_2$ + H which would otherwise convert fuel radicals directly to major products. The third is a constraint on alkylperoxides, APO ≡ CH${}_3$OOH + CH${}_3$OO + CH${}_2$OOH, imposed by slow reactions which convert APO to hydroperoxides coupled with fast reactions which equilibrate the species comprising APO. The fourth is a constraint on alcohol plus formaldehyde, ALCD ≡ CH${}_3$OH + CH${}_3$O + CH${}_2$OH + CH${}_2$O imposed by relatively slow reactions which generate/remove ALCD coupled with fast reactions which equilibrate the species comprising ALCD.

The following two combustion examples show typical RCCE results under the constraints discussed above, in comparison with the corresponding DKM calculations.

Given the fact that the internal dynamics of the system and the controlling constraints depend, in general, upon the initial distance of the system from its chemical equilibrium state, we consider two extreme cases below: close-to and far-from equilibrium systems.

5.1. Example 1. Close-to-Equilibrium Oxygen-Methane Combustion

The constraints defined in the previous section are particularly useful for modeling systems driven away from an initial equilibrium by fast expansions and compressions. They are useful both to increase the speed and efficiency of calculations and to provide insight about the most important reaction mechanisms. An example of this is shown in Figure 1 where RCCE calculations of the CO mole fraction during the expansion of combustion products in an internal combustion engine are compared with a detailed (DKM) calculation [35]. RCCE calculations are carried out using rate equations for the constraint-potentials.

Figure 1. (Color online) CO mole fractions as a function crank angle during the expansion stroke of an internal combustion engine.

Figure 1. (Color online) CO mole fractions as a function crank angle during the expansion stroke of an internal combustion engine.

RCCE calculations involving only energy and volume constraints are shown by the curve labeled LTE, i.e., the shifting equilibrium case. The DKM calculation includes 29 species and 132 reactions abstracted from the GRI-3 mechanism [80] and the compilations in [81,82]. It can be seen that the slow rate of expansion during the early stages of the power stroke results in local adjustment of composition to LTE predictions. However, as the rate of expansion increases, the LTE prediction departs significantly from the DKM prediction, due to the finite rates of some of the kinetic mechanisms that determine the relaxation process. The LTE results at the end of the expansion stroke are lower than the DKM results by about an order of magnitude. Also shown in Figure 1 is the choice-of-constraints dependence study of the RCCE predictions. Consistently with the Le Ch$\widehat{\mathrm{a}}$telier-Braun principle, the shift of the internal dynamics towards the exothermic direction due to the cooling effect of expansion makes the recombination reactions an essential part of the energy-restoration dynamics. Just constraining the total number of moles of gas, M, shows a significant improvement with respect to the LTE predictions. Adding constraints one at a time shows that only a subset of structural constraints, namely M, FV, and FO, are sufficient to capture the nonequilibrium dynamics during expansion. Again, we note that the constraint-potentials obtained from the RCCE calculations can be used to obtain mole fractions not only for the 29 species in the detailed model but also for any species involving the same elements and for which the thermodynamic properties are known.

5.2. Example 2. Far-from-Equilibrium Oxygen-Methane Combustion

For far-from-equilibrium systems, it is necessary to use the fuel as a constraint to fix the initial state of the system. The fuel constraint together with the elemental constraints are sufficient to guarantee that the system will go to the correct final equilibrium state. However, to obtain the correct time evolution of the system, additional constraints may be needed and the accuracy of the results depends on the careful selection of these added constraints.

This case is illustrated in Figure 2 which shows the predicted temperature profiles when constraints are added one at a time for a homogeneous stoichiometric mixtures of CH${}_4$/O${}_2$ in a constant volume adiabatic chamber. The initial temperature is 900 K and the initial pressure 10 atm. The DKM calculations include 132 reactions and 29 species, and are again based on the GRI-Mech 3.0 and the other compilations previously cited. The minimum number of reactions needed in the RCCE model, however, depends on the number and the type of constraints. Nevertheless, since all such calculations involve elemental constraints, the constrained equilibrium concentration of all 29 species can be obtained. The improvement in the prediction of the temperature profile can be seen to steadily increase over the entire range of conditions studied as additional constraints are included in the model.

In general, the accuracy of RCCE calculations improves as the number of constraints increases. In the limit where the number of constraints approaches the number of species in the DKM, the RCCE and DKM calculations become almost equivalent, but RCCE allows to evaluate concentrations also of species not included in DKM.

Although strongly supported by its general thermodynamic consistency, the experience and intuition dependent “art” of choosing the most appropriate constraints for an RCCE calculation may not be as effective when experimental data and reliability of the DKM are scarce. In such cases, we suggest that systematic studies of the DKM based on the tools of dynamical systems theory, such as the eigenvalue analysis of the local Jacobians [83] and the construction of the slow one-dimensional invariant manifold (as done in the SIM, ILDM, CSP, MIM, G and other schemes cited in the Introduction), could provide useful hints for the identification of physically meaningful structural constraining functions that maintain their rate-controlling capabilities for a broad class of problems, independently (or for a broad range) of initial thermodynamic conditions.

Figure 2. (Color online) Constraint-dependence study of RCCE predictions of the temperature profile during the oxidation of stoichiometric methane-oxygen mixtures at $T_i=900$ K and $p_i=10$ atm in a constant volume adiabatic chamber.

Figure 2. (Color online) Constraint-dependence study of RCCE predictions of the temperature profile during the oxidation of stoichiometric methane-oxygen mixtures at $T_i=900$ K and $p_i=10$ atm in a constant volume adiabatic chamber.

6. Conclusions

RCCE is a powerful method for simplifying the description of non-equilibrium states and their tendency to relax towards equilibrium. In this paper, we present its generalization in a generic mathematical context where the constraints can be nonlinear functions of the state variables. In doing so, we review the philosophy of the method, we discuss in a careful step-by-step derivation the various assumptions that underlie its application in chemical kinetics, and by showing the results of two typical combustion examples we discuss the basic structural constraints which successfully provide accurate model reduction for two different reacting systems.

We show that in the thermodynamics and statistical mechanics contexts, RCCE is fully compatible with the Second Law. In fact, it is a logical extension of conventional equilibrium thermodynamics to nonequilibrium systems. The method has the potential to greatly simplify the treatment of complex reacting systems for which only incomplete data are available about the reaction rates of the underlying detailed scheme, which often must be modeled by truncating the species list or introducing ad hoc global reactions. It can also be used to build combustion models for complex chemical systems for which standard methods to compute the rate constants are not accurate enough to justify individual-reaction modeling.

Among the relevant features of the RCCE method in the chemical kinetics context are the following:

• By focusing attention on the rate-limiting reactions, the method leads to a better understanding of the important mechanisms and reaction paths involved in the evolution of a complex system.
• Since the number of constraints required for RCCE calculations is much smaller than the number of species required for detailed calculations, the total number of differential plus algebraic equations to be solved is much smaller than for other models, so that the model is suitable for coupling with transport equations and turbulence models for turbulent flow combustion calculations.
• Since only reactions that change the constraints are needed, the number of reaction rates and the level of effort required to tabulate and input them are greatly reduced.
• RCCE combustion calculations can be initiated with only the constraints on the state variables and the fuel, and then systematically improved by adding constraining functionals of the state variables, one at a time, until an acceptable level of accuracy is reached. Excellent agreement with experimental data and validated detailed kinetic model calculations is often reached with a number of constraints well below the number of species in a detailed model.
• The entropy production rate in RCCE calculations is always non-negative. Moreover, the correct final equilibrium state for the fixed atomic-nuclei constraints is always reached, regardless of the number of constraints employed.

The present generalized mathematical formulation extends the applicability of the RCCE method to a wider class of relaxation problems and nonequilibrium dynamical systems, possibly even beyond the framework of chemical kinetics, provided they admit a Lyapunov functional of the state variables which cannot be increased by the system’s internal dynamics.