- Email: [email protected]
A hierarchical approach for fast and accurate equilibrium calculation
Computer Physics Communications 163 (2004) 177–190 www.elsevier.com/locate/cpc
A hierarchical approach for fast and accurate equilibrium calculation G. Colonna a,∗ , A. D’Angola b a CNR-IMIP, Bari section and Dipartimento di Chimica, Universitàdi Bari, via Orabona, 4, 70126 Bari, Italy b Dipartimento di Ingegneria e Fisica dell’Ambiente, Universitàdella Basilicata, C.da Macchia Romana, 85100 Potenza, Italy
Received 15 February 2004; received in revised form 2 August 2004; accepted 24 August 2004 Available online 1 October 2004
Abstract In the present paper an improved stoichiometric algorithm based on the successive single reaction equilibrium has been described, avoiding the solution of large systems. The algorithm calculates the equilibrium composition adding automatically suitable new reactions to improve the convergence. The algorithm has been tested for plasma composition calculation for a pure nitrogen and air plasma mixtures. The algorithm proposed always converges and it is significantly faster than other hierarchical algorithms. 2004 Elsevier B.V. All rights reserved. PACS: 82.60.Hc; 52.25.Kn Keywords: Homogeneous gas phase equilibrium; Chemical equilibrium; Method of calculation; Plasma composition; Stoichiometric algorithm
1. Introduction The numerical calculation of chemical equilibrium started in the fifties but is still widely studied. The algorithms developed can be classified in different ways [1]: stoichiometric or non-stoichiometric approaches, Gibbs free energy minimization [2] or mass action low solution, simultaneous solution of all the equi* Corresponding author. Tel.: +390805442082, fax: +390805442024. E-mail addresses: [email protected] (G. Colonna), [email protected] (A. D’Angola).
0010-4655/$ – see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.08.004
librium equations or solving one reaction at time hierarchically. This last method has been almost forgotten and the only approach survived is the solution of non-linear systems of equations based on Newton– Raphson (NR) or similar algorithms. What is strange is that hierarchical approach is more straightforward and easier such that it is used as hand calculation method for simple multireaction equilibrium. In fact one of the first algorithms, proposed by Villars [3, 4], uses the hierarchical approach. After these works, Cruise [5] used an approach similar to Villars’ [3,4], while Meintjes et al. [6] combined Villars’ approach with the NR method. Nowadays the method was com-
178
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
pletely disappeared from the literature. Villars himself modified his approach to solve simultaneously all the reactions. The advantage of solving one reaction at time is that “it does not require the inversion of large matrices” [5], speeding up the convergence with respect to NR by a factor 5. Villars [3] affirms that his algorithm always converge, but the speed of convergence depends strongly on the initial conditions. Therefore speeding up the convergence depends on the capacity of the researcher in finding the right initial conditions. He proposed also a way to improve the algorithm: adding some ad hoc reactions can speed up the convergence. The main drawbacks of the Villars’ method is that in cases when the convergence is slow it is not possible to speed up the convergence automatically, but the researcher should find on his experience how to improve the convergence. On the other hand matrix handling is straightforward for modern computer, and the solution of nonlinear systems and minimization techniques has been considerably improved in the last years because of the interest in many fields. The synergy between Villars’ drawbacks and matrix algorithm improvements, makes the hierarchical approach to disappear. In this work we want to propose a modified Villars method which shows a very fast convergence. Even if the general aspects of the method here presented are quite similar to Villars’, the drawbacks have been eliminated. First of all we do not need to consider an “optimum basis” as in [5]. We have a different criterion to choose the order of reactions that gives a better convergence. Moreover we have an automatic algorithm to introduce new reactions that speed up the convergence. The resulting algorithm is a completely automatized efficient version of the hierarchical approach, independently on the system studied.
2. General aspects In this paragraph we will discuss a numerical approach to calculate the equilibrium composition of homogeneous systems. In particular we will describe gas phase (plasma) equilibrium. Anyway, the method here described is based on the solution on algebraic equations and therefore it can be easily extended to solution of heterogeneous reactions. There are many different
aspects that must be considered to use this new approach, from the formalism to numerical precision. Any of this part is critical to have an efficient equilibrium solver. 2.1. Formalism To define a “chemical reacting system” we need to determine the set of chemical species that can be present in the system, the chemical reactions with their equilibrium constants and the initial composition. Here for initial composition we intend the mass ratio of elementary components (atoms or radicals), being this ratio constant for a reacting system. For gas phase systems there are two other quantities that must be fixed, one that determine the energy content (i.e. temperature, enthalpy or entropy), and one that gives the mass content (pressure or mass density). These quantities give two constraint equations to the reaction equilibrium. The chemical reactions will be written in the general form nrj
np
(j ) Cr j,i Ri
⇔
i=1
j
(j )
Cpj,i Pi ,
(1)
i=1
where nrj and npj are the number of reactants and (j )
(j )
products for the j th reaction, Ri and Pi are the reacting species and Crj,i and Cpj,i are the stoichiometric coefficients, that are all positive. An equivalent way of writing the reactions, which is the standard for traditional solvers, is nj
Cj,i Xi = 0,
(2)
i=1
where X’s are the set of all species and the C’s are the stoichiometric coefficients that can be positive (for reactants) and/or negative (for products). In our approach, reactants and products must be treated differently and therefore expression (1) is more practical, even if there is no real need to use representation in Eq. (1) instead of the one in Eq. (2). Following the scheme in Eq. (1), given the number of reactants and products, we store for each reaction the index of reactants and of products with their positive stoichiometric coefficients separately. Given a set of species concentrations {[Xi ]}, it is possible to define a quantity for each reaction called
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
reaction quotient (Πj ) npj (j ) Cp [Pi ] j,i Πpj = . Πj = i=1 nrj (j ) Crj,i Πrj i=1 [Ri ]
(3)
The unit of measure of the species concentration [Xi ] can be the molarity, the molar fraction or the partial pressure. The mass action low (equilibrium equation) can be written as the (eq)
Kj = Πj
,
(4)
where “Kj ” is the appropriate equilibrium constant of the j th reaction and the index (eq) refers to quantities calculated at equilibrium conditions. It must be noted that independently on the quantities used to express the species concentration (partial pressures, moles per unit volume and so on), the formalism and the equations does not change if one uses the right equilibrium constant (Kp , Kc , etc.). In this paper we prefer to use partial pressures, but the equilibrium routines can work with any kind of composition quantities.
a concentration values set (labeled with a “0”), for a single reaction, the system has only one degree of freedom to conserve the ratio between the atomic component masses. Therefore only one real variable x, called extent of reaction, is needed to define the state of the one reaction system (j ) (j ) Ri = Ri 0 − Cr j,i · x, (j ) (j ) Pi = Pi 0 − Cpj,i · x. As a consequence Πj , Πrj , Πpj are function of x and in particular Πrj , Πpj are polynomial of x while Πj is a rational function (see Eq. (3)). Therefore, the following function Fj (x) = Kj · Πrj − Πpj
(7)
i=1
It must be noted that the function F has the same properties as the reaction position ϑ.
Let’s consider a system with only one reaction, even if, for sake of generality, in the following we will retain the reaction index j . Eq. (4) holds only if the concentrations of the species have reached the equilibrium values. Out of equilibrium Πj can be lower or higher than Kj , corresponding respectively to reactant or product concentrations higher than the equilibrium values. On these bases, given any set of concentration values, it is possible to define a quantity, that we call reaction position (with respect to equilibrium) and indicated as ϑ, which has the following properties: (a) ϑ > 0 if the reaction is shifted toward reactants; (b) ϑ < 0 if the reaction is shifted toward products; (c) ϑ = 0 if the reaction is at equilibrium. There are different possible definitions of the quantity ϑ that we will discuss later. As a matter of fact, all these definitions must be equivalent and, as a consequence of Eqs. (3), (4), they must depend on the quantity Kj − Πj .
(6)
is a polynomial of degree nrj npj deg(Fj ) = max Crj,i , Cpj,i . i=1
2.2. Reaction position and distance
179
(5)
Moreover, it is possible to define a distance from equilibrium δ simply as the absolute value of ϑ. Given
2.3. Single reaction equilibrium To calculate the equilibrium composition we apply Eq. (4) to the function (6) obtaining that the x value to calculate the equilibrium composition (xeq ) from Eq. (4) is the solution of the equation Fj (xeq ) = 0.
(8)
The main problem in calculating equilibrium composition is that Eq. (8) has more than one solution. Therefore it is necessary to find a criterion that automatically choose the physical solution between all the possible ones. The approach used consists in finding an interval in which only the physical solution of Eq. (8) is included. The first extreme of such interval is given by the value x0 = 0 (initial conditions). To determine the other extreme we use the properties of ϑ listed above: its sign gives the reaction direction. If ϑ > 0 the reactant concentrations are higher than equilibrium and therefore the other extreme is obtained for a positive value of x. The limiting value can be estimated supposing that the reaction is complete, so that at least one reactant has vanishing concentration, which can
180
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
be expressed by the relation (j ) [Ri ]0 xmax = min . Cr j,i
(9)
Similarly if ϑ < 0, the other extreme of x is negative and is evaluated as (j ) [Pi ]0 . xmin = − min (10) Cpj,i The existence of a solution in such interval is assured by the fact that the function F has opposite sign in the two extreme points. The unicity of the solution in such interval can be easily demonstrated mathematically. Another problem in finding the equilibrium composition is to solve Eq. (8) efficiently and with a good accuracy considering that the polynomial coefficients can range between many orders of magnitude. In the general case it must be calculated by iterative algorithm (NR, secant, etc.). The NR algorithm is not always stable so that one must switch to the other methods when NR does not converge or converge to a non physical solution. It is also difficult to decide which is the best value of convergence parameters nor even the convergence criterion due to the problem stiffness. Moreover the precision required in the solution of the equation must be very high because the reaction distance can be sensitive to the solution. For this reason a particular care must be taken to decide the convergence of root finding routines. We have developed an algorithm based on the Newton, secant and midpoint methods, coupled to improve the convergence. Also the secant method shows some problems due to truncation errors. If the solution is very close to one of the extreme, the convergence can be very slow, but it is possible to improve the convergence simply considering the new solution as a weighted mean of the interval extremes giving the higher weight to the one closer to the solution. This method is not always better than the pure secant, but it always keeps the convergence to a reasonable speed. When the secant method fails to decrease the interval amplitude, we uses the midpoint rule, which always converges. The criterion for convergence considers that the value of the function in x must be lower that a tolerance value (Tols). This criterion fails when truncation errors are important. Therefore an additional criterion can be
added: if the interval amplitude is lower than a given value, comparable to the round off, the convergence has been reached. Differently from [3,5] where all the reactions are expressed starting from a species basis, in our approach we prefer to include, as much as possible, reactions with no more than two reactants and/or products, which produces polynomial of second order degree for Eq. (8), that can be resolved directly. It is important to note that, the stoichiometry of the system is always conserved. Moreover the solution of the equation never gives negative concentrations due to the limiting conditions (see Eqs. (9), (10)). Non physical solutions have never been met, but round off errors can give a null value of x. This problem has been avoided calculating approximate solution for x. For a second degree equation ax 2 + bx + c = 0, if a · c b2 , below the round off error, the direct solutions of the equation are x1 = 0 and x2 = −b/a, as in the case with c = 0. This problem can be easily solved remembering that x1 · x2 = c/a. Because x2 can be considered precise to the last digit, x1 can be approximated as x1 = c/a/x2 = −c/b which correspond to the solution obtained neglecting the terms with the power of x higher than one. This result can be easily extended to polynomial of degree higher than two. 2.4. Multi-reaction systems The calculation of equilibrium for multi-reaction systems could start from the equilibrium of a single reaction. The method here proposed consists in calculating the equilibrium concentrations solving hierarchically one-reaction equilibrium problem at a time, until the convergence is reached. Probably, the algorithm converges independently on the sequence of reactions chosen in the calculation [3], but it can be optimized choosing an appropriate sequence. In most cases, some reactions can be neglected, because they are almost ineffective for the equilibrium composition. At this purpose, the reaction distance from the equilibrium can be a suitable parameter to choose the reaction order: higher the distance, further the reaction from the equilibrium. The algorithm for the solution of the equilibrium problem for multireaction systems can be organized in the following way: starting from an initial value of concentration, we calculate the reaction position (ϑ) and distance (δ) for all the reactions. Than
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
we calculate the equilibrium composition for the reaction with the maximum distance. We iterate the previous procedure starting from the last calculated concentrations values. The solution of the single reaction equilibrium must be calculated with a good accuracy, otherwise, due to the large differences in the equilibrium constants, the reaction just considered could be solved again uselessly, reducing efficiency and accuracy of the overall equilibrium calculation. As stated before, there is more than one possible reaction position. We have chosen two different definitions. They behave differently and present advantages and drawbacks. The first definition (ϑ1 ) is ϑ1j = Fj (0) and the second (ϑ2 ) Kj · Πrj − 1 Kj · Πrj < Πpj , Πpj Πpj ϑ2j = Kj · Πrj > Πpj , 1− Kj · Πrj 0 Kj · Πrj = Πpj .
(11)
(12)
Both the definitions have the advantage to be continuous functions of the extent of reaction and finite real numbers even if some reactants or products have zero concentration. Anyway they behave differently: ϑ1 is not limited, while ϑ2 is between −1 (when at least one reactant has vanishing concentration) and 1 (when at least one product has vanishing concentration). For this reason the two positions are useful in different situations. The following example will clarify the problem: consider a diatomic molecule A2 that can be dissociated (a: A2 → 2A) and ionized (b: − A2 → A+ 2 + e ). Starting with the condition [A2 ]0 = 1 + and [A2 ]0 = [e− ] = [A] = 0, we have that the first reaction positions are ϑ1a = Ka and ϑ1b = Kb and the second are ϑ2a = ϑ2b = 1. Therefore ϑ2 cannot be used to discriminate the two reactions and therefore the choice can be made only on random bases. Suppose, for example, that the A+ 2 concentration is very small compared to the atom concentration and that the first reaction chosen is the ionization process. What happens is that after the ionization, the dissociation reaction looses some molecules, changing the ionization equilibrium toward the recombination weakly, so that the convergence is very slow. On the contrary ϑ1 chooses the reaction with higher equilibrium constant, which is the reaction furthest from the equilib-
181
rium. The drawbacks of using ϑ1 are: the round off errors of reactions with high equilibrium constants can be higher than the value of ϑ1 for reactions very far from equilibrium with small equilibrium constants and therefore such reactions will never be considered. As a consequence, the solution obtained gives the right concentration only for majority species, while secondary species present large relative errors. To take advantage of both distances, we build up two routines, with the same algorithm but using different distances. Firstly we call the routine that uses ϑ1 till the convergence is reached. At this step, reactions with small equilibrium constant have not been considered. Therefore we call the routine with the position ϑ2 . This distance does not discriminate the reactions by the equilibrium constant, therefore the reactions neglected in the first step, because of the small equilibrium constant, are now calculated determining the species with low concentrations. Obviously, all the reactions are included in the calculation for both distances; so the algorithm decides the reactions to consider. The calculation for each distance should terminate when all the reaction distances are below the tolerance required. While for ϑ2 it is trivial to satisfy this condition, for ϑ1 some troubles arise, because it is difficult to define an adequate tolerance value for all the reactions due to the large differences in the equilibrium constants. At this purpose an additional convergence criterion can be adopted based on the relative concentration error. Moreover, it is possible to decide if ϑ1 has reached the precision limit by considering ϑ2 which is a relative evaluation. This means that the in the first step, when the reaction position is ϑ1 , we calculate also ϑ2 ; if it is smaller than a given tolerance, ϑ1 is of the order of the round off error and then must be considered null. 2.5. Addition of new reactions Even if the double distance approach speeds up the calculation of equilibrium composition, some problems cannot be avoided: for mixtures with many species and many reactions, some numerical traps can be found and as a consequence the convergence slows down. This is the case in which two species concentration (A, B in Eq. (13)) are coupled by successive reactions passing by an intermediate species (X) with
182
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
much lower concentration Ra : Rb :
A + ··· ⇒ X + ···, X + ··· ⇒ B + ···.
(13)
In these conditions each reaction produces very small changes with respect to the absolute variations of A and B species, and therefore the two reactions in Eq. (13) are repeated many times one after the other till the global transfer is obtained. As an example, the dissociation/ionization described in the previous section can be considered, where A stands for the atom, B for the molecular ion and X for the molecule. The solution of the problems can be obtained considering a new reaction as the combination of the two written in Eq. (13) obtained eliminating the common element X as in Ra,b :
A + ··· ⇒ B + ···.
(14)
Villars [3] already pointed out this problem, suggesting as solution the addition of a new reaction. This criterion has not been generalized and is based on the experience of the researcher. For this reason the Villars approach has been abandoned. In this paper we propose an automatic method to add as many reactions as needed with the result to improve considerably the convergence. It must be noted that, differently by matrix methods, hierarchical approach can use any number of reactions. To choose the reactions to combine we count the number of time each reaction has been solved after a given number of iteration. If one of these reactions is higher than a threshold, it is combined with the reaction with the higher number of repetition which matches the requirements in Eq. (13). In practice, adding a new reaction is finding a shortcut to the path that the system must follow to reach the equilibrium. The complete algorithm for calculating the equilibrium with the addition of new reactions works in this way (see the flow chart in Fig. 1). We have two nested loops. The external loop (a “while” loop) checks only if the convergence is reached, the internal loop (a “for” loop) repeats the calculation of single reaction equilibrium a given number of iterations (Nint ). Within the internal loop we count the number of time each reaction has been chosen. If a given number of repetition (Nrep ) is reached by a reaction, we check between the other reactions with the higher repetition number if they fulfill the requirement in Eq. (13) so that they
can be combined as described in Eq. (14). If two reactions are combined in a new one, we solve this new reaction and we reset the reaction counters of the combined reactions. When the internal loop is completed, we check if the convergence requirements are reached and we reset the counters for all reactions. The values chosen for Nint and Nrep weakly affect the convergence behaviors, but they must clearly fulfill the requirement Nint /Nrep 2, otherwise the condition to build the new reaction will never be met. As stated before, we decide to build the new reaction when two reactions are repeated one after the other as described in the dissociation–ionization example. This is not the only situation can give problems. A more longer series of reactions can be looped together giving similar troubles. The probability of finding this situation decreases while the number of the reaction series increases. The kind of series we can get with the nested loop technique depend on the ratio Nint /Nrep . 2.6. Constraints A short comment about constraints is needed. Constant pressure and constant volume constraints, or any other referring to the number of molecules in the system, can be easily included in the code introducing an appropriate normalization condition after each equation has reached the equilibrium. For example, for an ideal gas at constant pressure P , the partial pressures are normalized by the factor P / i Pi , while the constant volume condition does not require any normalization. For non ideal gases such as plasma, the equation of state must be corrected according to the Debye–Huckel (DH) theory [11] 1 = Pid + PDH , P = Kb T n − (15) 24πλ3d where λd is the Debye length, n the total number density of particles and Kb the Boltzmann constant. In the frame of this theory, the correction to the pressure depends on the charged particle density. As for the number of molecule constraint, also for the DH plasma, it is sufficient to write an appropriate normalization condition. For example, for the constant pressure DH plasma we should normalize the partial pressures by a
factor (P − PDH )/ i Pi . The convergence of the DH corrections is faster than equilibrium, because it is important only if the ionization degree is larger than or
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
Fig. 1. Flow chart of the equilibrium calculation algorithm.
183
184
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
close to one and this algorithm determine the majority species concentration very fast and accurately. We must say that the DH theory modify the equilibrium constants as well. This is not a problem in the model because the convergence of the DH correction is very fast and not very sensitive to charge concentrations. We must point out that the DH equation of state is in the perturbation theory framework, which means that the deviation from the ideal gas is small. We do not apply the algorithm to equations of state with strong perturbations. The constant temperature constraint obviously gives no problem while the constant energy one introduces a further complication due to the dependence of the equilibrium constants on the temperature. This problem can be treated in different ways but we do not have a definitive answer on it. The temperature can be updated after all the equilibrium has been calculated or during the calculation of equilibrium after each reaction or a given number of reactions has been solved. We choose this second approach, but we did not test the convergence properties in one case or the other.
putting in evidence the improvement due to the added reactions. At this purpose we have applied the method to ionized gaseous mixtures in the temperature range 1000–30000 K. In particular we have considered a pure nitrogen mixture (see the first three reactions in Table 1 and Fig. 2) and the air mixture (see Table 1 and Fig. 3). The first system is very simple, therefore we can investigate the behavior of the algorithm step by step, being each step (iteration) the resolution of a single chemical reaction. It must be observed that the results obtained for this system do not include multiple ionized particles even if at high temperatures they
3. Results and discussion In this section we will discuss the results obtained by the method described above in different cases,
Fig. 2. Equilibrium composition of a nitrogen mixture as a function of temperature.
Table 1 Equilibrium constants of the reactions considered at different temperature. Species considered: N2 , N, N+2 , N+ , N++ , N+++ , N++++ , O2 , O+2 , − + ++ +++ , O++++ , NO, NO+ , e− O− 2 , O, O , O , O , O Reactions
log(Keq )
(16)
103 K
5 · 103 K
104 K
2 · 104 K
3 · 104 K
N2 2N N2 N2 + e− N N+ + e− N+ N++ + e− N++ N+++ + e− N+++ N++++ + e− O2 O+2 + e− O2 + e− O− 2 O2 2O + − OO +e O+ O++ + e− O++ O+++ + e− O+++ O++++ + e− N2 + O2 2NO NO NO+ + e− O + e− O−
−85.93590 −165.3994 −153.5874 −329.5341 −537.5887 −883.7819 −125.6238 −7.681201 −32.49372 −144.1756 −392.8815 −623.5494 −885.2708 −18.65947 −94.88367 3.713479
5.036874 −16.59715 −14.39963 −50.56606 −93.27721 −160.5291 −9.595570 −16.56471 15.62565 −13.81455 −62.25641 −109.3668 −162.8089 −1.346583 −4.713501 −15.26020
16.94815 3.624266 4.127968 −14.51822 −36.50893 −68.89310 5.960819 −19.04409 21.72691 3.714115 −19.72256 −43.89081 −71.27454 0.6411276 7.757746 −18.91194
24.26754 14.49052 13.53105 4.296274 −7.222843 −22.26790 14.84219 −21.37386 25.42207 13.35034 2.214520 −10.36335 −24.62963 1.521508 15.15033 −21.88424
29.62973 18.13983 15.70289 10.86431 3.104701 −6.460977 18.37340 −22.68414 28.97627 16.15614 9.825148 1.235014 −8.567586 1.581945 18.43401 −24.35636
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
185
Fig. 3. Equilibrium composition of air mixture as a function of temperature.
become important. In the second one we want to investigate the behavior of the algorithm for realistic large systems. In Table 1 we have reported also the equilibrium constants for each reaction at different temperatures showing the stiffness of the equilibrium problem. The calculation of the equilibrium constants is made starting from energy levels and partition functions. There are different aspects that need to be verified. In particular we consider the double distance approach and the addition of a new reaction. We consider four cases which differ by the reaction position considered, discussing the convergence: (1) the distance is ϑ2 . The first reaction solved is molecular ionization (reaction 2 in Table 1); (2) the distance is ϑ2 . The first reaction solved is the molecular dissociation (reaction 1 in Table 1); (3) the distance is ϑ1 . (4) combined solution of ϑ1 and ϑ2 .
The evolution of the molar fractions toward the convergence values with the number of reactions solved for pure N2 mixture has been reported in Figs. 4–8 at gas temperature T = 30000 K. As initial composition molecular nitrogen at p = 1 atm has been considered, which is very far from the equilibrium composition. The equilibrium composition has been calculated at constant volume. In Figs. 4a, b it is reported the gas composition as a function of the number of reactions solved in case 1. To reach the convergence (gray lines) all the species needs more than 10000 iterations. The N2 concentration (Fig. 4b) presents strong oscillations. This is exactly equivalent to the example of dissociation/ionization system discussed in Section 2.4. In this case the first reaction solved is the N2 ionization. As a result, the gas is made mainly by ionized molecules. Therefore the second reaction solved is the dissociation. These steps push the solution farther from the equilibrium, being N+ the dominant species. For
186
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
(a)
(b)
Fig. 4. Evolution of the molar fractions with the number of reactions solved for N2 system at T = 30000 K, case 1. The gray lines represent the equilibrium values at convergence.
(a)
(b)
Fig. 5. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 2. The gray lines represent the equilibrium values at convergence.
this reason the system goes back to the equilibrium changing slowly the ionized molecule and atomic nitrogen concentrations, alternating the solution of N2 ionization and dissociation. Each couple of reactions changes slowly the concentrations (N and N+2 ) be-
cause the intermediate species (N2 ) is negligible. After 10000 reactions solved the system is still very far from equilibrium and ionization and dissociation reactions are still solved alternatively making N2 concentration to oscillate on an interval of some orders of magnitude.
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
(a)
187
(b)
Fig. 6. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 3. The gray lines represent the equilibrium values at convergence.
The convergence changes if the first reaction solved is the dissociation (case 2), as can be observed in Figs. 5a, b. The majority species (N+ , N) reach the equilibrium values after few iterations (3 or 4), while minority species (N2 , N+2 ) need more than 1000 steps to converge. Also in this case there are strong oscillations of the N2 convergence series due to the alternate solution of molecule dissociation and ionization as in case 1. These results show the limit of ϑ2 : it is not appropriate to decide which is the most important reaction. If the right choice is made at the first step, the convergence can be speeded up of one order of magnitude. This problem disappears using ϑ1 (case 3, see Figs. 6a, b). In fact, at the beginning, the reaction farthest from the equilibrium is the one with highest equilibrium constant, in this case the dissociation. By using ϑ1 , the convergence values of majority species are reached after few iterations while minority species concentrations converge slowly (more than 1000 iterations), showing the same behaviors as case 2. This is a very important result, because, using ϑ1 , if only high concentration species are relevant, the convergence can be reached rapidly. In the last case (4, see Fig. 7) using ϑ1 the system needs few iterations to reach the convergence, but mi-
Fig. 7. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 4. The gray lines represent the equilibrium values at convergence.
nority species present wrong concentrations. On the other hand, using ϑ2 , also minority species concentrations are calculated with the required tolerance, even if they are very small. The oscillations in N2 concentration observed in all the cases discussed above put in evidence the presence
188
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
of some troubles in this algorithm. Adding new reactions introduces a shortcut to the solution, eliminating the oscillating behavior. The set of reactions used in the previous section is minimal, i.e. any further reaction is a linear combination of the minimal set. This approach is necessary for standard methods which needs to solve linear systems while non-minimal set leads to overdimensioned systems. In our algorithm only one equation at a time is solved, therefore there is no problem in using non-minimal sets. The problem is then to choose the set of reactions which leads to the equilibrium solution as fast as possible but this
Fig. 8. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 1 with the addition of new reactions. The gray lines represent the equilibrium values at convergence.
choice is not possible a priori, because it depends on the thermodynamic conditions. In Fig. 8 we have reported the results obtained in case 1, adding new reactions automatically. In about 20 reactions solved, the equilibrium compositions are reached for all the species. It must be noted that there are only few oscillations in N2 concentrations, differently from Fig. 4b, because with the aid of the added reactions the concentrations jump directly to the correct values. In Tables 2, 3 we have reported the number of reactions to reach the convergence for a required tolerance (10−13 ) in the four cases listed above with (+R) and without new reactions for different temperatures,
Fig. 9. Evolution of the absolute value of ϑ2 with the number of iterations for air mixture at 1 atm and T = 10000 K, in different cases.
Table 2 Number of reactions solved to reach the convergence in different cases (see Section 3) and for different temperatures with (+R) and without the addition of further reactions for the pure nitrogen mixture. The tolerance required is 10−13 Cases
T = 103 K
T = 5 · 103 K
T = 104 K
T = 2 · 104 K
T = 3 · 104 K
1 1+R
3 3
31 36
31 27
1396 35
46506 24
2 2+R
3 3
30 28
29 33
818 31
24066 19
3 3+R
1 1
37 41
27 35
564 15
9328 12
4 4+R
3 3
31 57
27 66
865 52
25851 22
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
189
Table 3 Number of reactions solved to reach the convergence in different cases (see Section 3) and for different temperatures with (+R) and without the addition of further reactions for the air mixture. The tolerance required is 10−13 Cases
T = 103 K
T = 5 · 103 K
T = 104 K
T = 2 · 104 K
T = 3 · 104 K
1 1+R
59 59
510 302
2768 294
22591 684
>200000 435
2 2+R
57 62
520 406
2592 309
14589 833
>200000 327
3 3+R
3 3
641 128
2832 414
1277 171
19038 107
4 4+R
51 51
467 172
2394 328
13926 393
>200000 337
where the plasma is in completely different conditions, respectively for the pure N2 and air mixtures. In Fig. 9 we have reported the evolution of the absolute value of ϑ2 with the number of iterations for air mixture at 1 atm and T = 10000 K, in the cases labeled 1, 4 and 4 + R. The decreasing of |ϑ2 | is exponential and moreover the addiction of new reactions speeds up the convergence rate. In this particular case, the case 1 + R is faster than 4 + R at low tolerance. This is not a general case as can be observed in Tables 2, 3. The main advantage of the automatic addition of new reactions can be observed when the equilibrium composition is far from the initial values, reducing the number of iterations of many orders of magnitude. It must be noted that independently on the case considered, the addition of new reactions retains the number of iteration in the same order of magnitude. It is worth noticing that using ϑ1 the algorithm converges faster. The drawback is that very low concentrations are miscalculated. Nevertheless in some applications, such as equilibrium fluid dynamics or transport properties calculations, this approach can be used because only majority species are important.
4. Conclusions In this paper we have presented an efficient stoichiometric, hierarchical method to calculate single phase equilibrium composition. The novelty of this method respect to previous hierarchical algorithms is the automatic addition of new reactions in order to improve the speed of convergence. We have studied two
gas mixtures in different conditions to validate the algorithm, analyzing the differences with traditional approaches. The results obtained show that the algorithm is very efficient in all the conditions, independently on the initial values. We have considered ideal gas approximation but the method has been also applied for non ideal plasma with the Debye–Huckel correction. The method is very suitable for non ideal plasma calculation where the correction terms depend on the free electron concentration. The automatic addition of new reactions makes the code very efficient in a wide range of conditions. We have used it to calculate air properties [12] in the temperature range 50–105 K and in the pressure range 10−3 –103 atm. In all these cases the code never fails differently from STANJAN [13] and NASA [14] whose convergence is related to the constraints used [2]. Moreover the convergence of the method is enforced by the fact that it never calculates at any iterations non physical solutions such as negative concentrations. This code is the first attempt after about 50 years [3] to propose a modern version of a hierarchical approach that is automatic and fast enough to be of practical use. Recent efforts on the chemical equilibrium calculation are devoted to multiphase systems [7–10]. All these methods are based on the simultaneously solution of the non linear equation mainly considering Gibbs energy minimization. In our opinion the method here presented can be easily extended to multiphase systems.
190
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
References [1] W.R. Smith, R.W. Missen, Chemical Reaction Equilibrium Analysis: Theory and Algorithms, John Wiley and Sons, New York, 1982. [2] P.S. Bishnuy, D. Hamirouney, M. Metghalchiy, J.C. Keckz, Constrained-equilibrium calculations for chemical systems subject to generalized linear constraints using the NASA and STANJAN equilibrium programs, Combust. Theory Modelling 1 (1997) 295. [3] D.S. Villars, in: G. Bahn, E. Zuckowsky (Eds.), Proceedings of the First Conference on the Kinetics Equilibria and Performance of High Temperature Systems, Western States Section of the Combustion Institute, Butterworths Scientific Publications, Washington DC, 1960. [4] D.S. Villars, A method of successive approximations for computing combustion equilibria on high speed digital computer, J. Phys. Chem. 60 (1959) 521. [5] D.R. Cruise, Notes on rapid computation of chemical equilibria, J. Phys. Chem. 68 (1964) 3797. [6] K. Meintjes, A.P. Morgan, Performance of algorithms for calculating the equilibrium composition of a mixture of gases, J. Comp. Phys. 60 (1985) 219.
[7] C.M. McDonald, C.A. Floudas, GLOPEQ: A new computational tool for the phase and chemical equilibrium problem, Comput. Chem. Engrg. 21 (1997) 1. [8] Y. Sofyan, A.J. Ghajar, K.A.M. Gasem, Multiphase equilibrium calculations using Gibbs minimization techniques, Ind. Eng. Chem. Res. 42 (2003) 3786. [9] A.V. Phoenix, R.A. Heidemann, A non-ideal multiphase chemical equilibrium algorithm, Fluid Phase Equil. 150–151 (1998) 255. [10] R.N. Mioshi, C.L. do Lago, An equilibrium simulator for multiphase equilibria based on the extent of reaction and Newton– Raphson method with globally convergent strategy (SEQEx2), Anal. Chim. Acta 334 (1996) 271. [11] L.D. Landau, Statistical Physics, Elsevier Science, 1979, Chapter 7. [12] M. Capitelli, G. Colonna, C. Gorse, A. D’Angola, Transport properties of high temperature air in local thermodynamic equilibrium, Eur. Phys. J. D 11 (2000) 279. [13] W.C. Reynolds, The element potential method for chemical equilibrium analysis: implementation in the interactive program STANJAN, Stanford University Report ME 270 HO no 7, 1986. [14] S. Gordon, B.J. McBride, Computer program for calculation of complex chemical equilibrium compositions and applications, NASA Reference Publication 1311 (1994).
A hierarchical approach for fast and accurate equilibrium calculation G. Colonna a,∗ , A. D’Angola b a CNR-IMIP, Bari section and Dipartimento di Chimica, Universitàdi Bari, via Orabona, 4, 70126 Bari, Italy b Dipartimento di Ingegneria e Fisica dell’Ambiente, Universitàdella Basilicata, C.da Macchia Romana, 85100 Potenza, Italy
Received 15 February 2004; received in revised form 2 August 2004; accepted 24 August 2004 Available online 1 October 2004
Abstract In the present paper an improved stoichiometric algorithm based on the successive single reaction equilibrium has been described, avoiding the solution of large systems. The algorithm calculates the equilibrium composition adding automatically suitable new reactions to improve the convergence. The algorithm has been tested for plasma composition calculation for a pure nitrogen and air plasma mixtures. The algorithm proposed always converges and it is significantly faster than other hierarchical algorithms. 2004 Elsevier B.V. All rights reserved. PACS: 82.60.Hc; 52.25.Kn Keywords: Homogeneous gas phase equilibrium; Chemical equilibrium; Method of calculation; Plasma composition; Stoichiometric algorithm
1. Introduction The numerical calculation of chemical equilibrium started in the fifties but is still widely studied. The algorithms developed can be classified in different ways [1]: stoichiometric or non-stoichiometric approaches, Gibbs free energy minimization [2] or mass action low solution, simultaneous solution of all the equi* Corresponding author. Tel.: +390805442082, fax: +390805442024. E-mail addresses: [email protected] (G. Colonna), [email protected] (A. D’Angola).
0010-4655/$ – see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cpc.2004.08.004
librium equations or solving one reaction at time hierarchically. This last method has been almost forgotten and the only approach survived is the solution of non-linear systems of equations based on Newton– Raphson (NR) or similar algorithms. What is strange is that hierarchical approach is more straightforward and easier such that it is used as hand calculation method for simple multireaction equilibrium. In fact one of the first algorithms, proposed by Villars [3, 4], uses the hierarchical approach. After these works, Cruise [5] used an approach similar to Villars’ [3,4], while Meintjes et al. [6] combined Villars’ approach with the NR method. Nowadays the method was com-
178
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
pletely disappeared from the literature. Villars himself modified his approach to solve simultaneously all the reactions. The advantage of solving one reaction at time is that “it does not require the inversion of large matrices” [5], speeding up the convergence with respect to NR by a factor 5. Villars [3] affirms that his algorithm always converge, but the speed of convergence depends strongly on the initial conditions. Therefore speeding up the convergence depends on the capacity of the researcher in finding the right initial conditions. He proposed also a way to improve the algorithm: adding some ad hoc reactions can speed up the convergence. The main drawbacks of the Villars’ method is that in cases when the convergence is slow it is not possible to speed up the convergence automatically, but the researcher should find on his experience how to improve the convergence. On the other hand matrix handling is straightforward for modern computer, and the solution of nonlinear systems and minimization techniques has been considerably improved in the last years because of the interest in many fields. The synergy between Villars’ drawbacks and matrix algorithm improvements, makes the hierarchical approach to disappear. In this work we want to propose a modified Villars method which shows a very fast convergence. Even if the general aspects of the method here presented are quite similar to Villars’, the drawbacks have been eliminated. First of all we do not need to consider an “optimum basis” as in [5]. We have a different criterion to choose the order of reactions that gives a better convergence. Moreover we have an automatic algorithm to introduce new reactions that speed up the convergence. The resulting algorithm is a completely automatized efficient version of the hierarchical approach, independently on the system studied.
2. General aspects In this paragraph we will discuss a numerical approach to calculate the equilibrium composition of homogeneous systems. In particular we will describe gas phase (plasma) equilibrium. Anyway, the method here described is based on the solution on algebraic equations and therefore it can be easily extended to solution of heterogeneous reactions. There are many different
aspects that must be considered to use this new approach, from the formalism to numerical precision. Any of this part is critical to have an efficient equilibrium solver. 2.1. Formalism To define a “chemical reacting system” we need to determine the set of chemical species that can be present in the system, the chemical reactions with their equilibrium constants and the initial composition. Here for initial composition we intend the mass ratio of elementary components (atoms or radicals), being this ratio constant for a reacting system. For gas phase systems there are two other quantities that must be fixed, one that determine the energy content (i.e. temperature, enthalpy or entropy), and one that gives the mass content (pressure or mass density). These quantities give two constraint equations to the reaction equilibrium. The chemical reactions will be written in the general form nrj
np
(j ) Cr j,i Ri
⇔
i=1
j
(j )
Cpj,i Pi ,
(1)
i=1
where nrj and npj are the number of reactants and (j )
(j )
products for the j th reaction, Ri and Pi are the reacting species and Crj,i and Cpj,i are the stoichiometric coefficients, that are all positive. An equivalent way of writing the reactions, which is the standard for traditional solvers, is nj
Cj,i Xi = 0,
(2)
i=1
where X’s are the set of all species and the C’s are the stoichiometric coefficients that can be positive (for reactants) and/or negative (for products). In our approach, reactants and products must be treated differently and therefore expression (1) is more practical, even if there is no real need to use representation in Eq. (1) instead of the one in Eq. (2). Following the scheme in Eq. (1), given the number of reactants and products, we store for each reaction the index of reactants and of products with their positive stoichiometric coefficients separately. Given a set of species concentrations {[Xi ]}, it is possible to define a quantity for each reaction called
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
reaction quotient (Πj ) npj (j ) Cp [Pi ] j,i Πpj = . Πj = i=1 nrj (j ) Crj,i Πrj i=1 [Ri ]
(3)
The unit of measure of the species concentration [Xi ] can be the molarity, the molar fraction or the partial pressure. The mass action low (equilibrium equation) can be written as the (eq)
Kj = Πj
,
(4)
where “Kj ” is the appropriate equilibrium constant of the j th reaction and the index (eq) refers to quantities calculated at equilibrium conditions. It must be noted that independently on the quantities used to express the species concentration (partial pressures, moles per unit volume and so on), the formalism and the equations does not change if one uses the right equilibrium constant (Kp , Kc , etc.). In this paper we prefer to use partial pressures, but the equilibrium routines can work with any kind of composition quantities.
a concentration values set (labeled with a “0”), for a single reaction, the system has only one degree of freedom to conserve the ratio between the atomic component masses. Therefore only one real variable x, called extent of reaction, is needed to define the state of the one reaction system (j ) (j ) Ri = Ri 0 − Cr j,i · x, (j ) (j ) Pi = Pi 0 − Cpj,i · x. As a consequence Πj , Πrj , Πpj are function of x and in particular Πrj , Πpj are polynomial of x while Πj is a rational function (see Eq. (3)). Therefore, the following function Fj (x) = Kj · Πrj − Πpj
(7)
i=1
It must be noted that the function F has the same properties as the reaction position ϑ.
Let’s consider a system with only one reaction, even if, for sake of generality, in the following we will retain the reaction index j . Eq. (4) holds only if the concentrations of the species have reached the equilibrium values. Out of equilibrium Πj can be lower or higher than Kj , corresponding respectively to reactant or product concentrations higher than the equilibrium values. On these bases, given any set of concentration values, it is possible to define a quantity, that we call reaction position (with respect to equilibrium) and indicated as ϑ, which has the following properties: (a) ϑ > 0 if the reaction is shifted toward reactants; (b) ϑ < 0 if the reaction is shifted toward products; (c) ϑ = 0 if the reaction is at equilibrium. There are different possible definitions of the quantity ϑ that we will discuss later. As a matter of fact, all these definitions must be equivalent and, as a consequence of Eqs. (3), (4), they must depend on the quantity Kj − Πj .
(6)
is a polynomial of degree nrj npj deg(Fj ) = max Crj,i , Cpj,i . i=1
2.2. Reaction position and distance
179
(5)
Moreover, it is possible to define a distance from equilibrium δ simply as the absolute value of ϑ. Given
2.3. Single reaction equilibrium To calculate the equilibrium composition we apply Eq. (4) to the function (6) obtaining that the x value to calculate the equilibrium composition (xeq ) from Eq. (4) is the solution of the equation Fj (xeq ) = 0.
(8)
The main problem in calculating equilibrium composition is that Eq. (8) has more than one solution. Therefore it is necessary to find a criterion that automatically choose the physical solution between all the possible ones. The approach used consists in finding an interval in which only the physical solution of Eq. (8) is included. The first extreme of such interval is given by the value x0 = 0 (initial conditions). To determine the other extreme we use the properties of ϑ listed above: its sign gives the reaction direction. If ϑ > 0 the reactant concentrations are higher than equilibrium and therefore the other extreme is obtained for a positive value of x. The limiting value can be estimated supposing that the reaction is complete, so that at least one reactant has vanishing concentration, which can
180
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
be expressed by the relation (j ) [Ri ]0 xmax = min . Cr j,i
(9)
Similarly if ϑ < 0, the other extreme of x is negative and is evaluated as (j ) [Pi ]0 . xmin = − min (10) Cpj,i The existence of a solution in such interval is assured by the fact that the function F has opposite sign in the two extreme points. The unicity of the solution in such interval can be easily demonstrated mathematically. Another problem in finding the equilibrium composition is to solve Eq. (8) efficiently and with a good accuracy considering that the polynomial coefficients can range between many orders of magnitude. In the general case it must be calculated by iterative algorithm (NR, secant, etc.). The NR algorithm is not always stable so that one must switch to the other methods when NR does not converge or converge to a non physical solution. It is also difficult to decide which is the best value of convergence parameters nor even the convergence criterion due to the problem stiffness. Moreover the precision required in the solution of the equation must be very high because the reaction distance can be sensitive to the solution. For this reason a particular care must be taken to decide the convergence of root finding routines. We have developed an algorithm based on the Newton, secant and midpoint methods, coupled to improve the convergence. Also the secant method shows some problems due to truncation errors. If the solution is very close to one of the extreme, the convergence can be very slow, but it is possible to improve the convergence simply considering the new solution as a weighted mean of the interval extremes giving the higher weight to the one closer to the solution. This method is not always better than the pure secant, but it always keeps the convergence to a reasonable speed. When the secant method fails to decrease the interval amplitude, we uses the midpoint rule, which always converges. The criterion for convergence considers that the value of the function in x must be lower that a tolerance value (Tols). This criterion fails when truncation errors are important. Therefore an additional criterion can be
added: if the interval amplitude is lower than a given value, comparable to the round off, the convergence has been reached. Differently from [3,5] where all the reactions are expressed starting from a species basis, in our approach we prefer to include, as much as possible, reactions with no more than two reactants and/or products, which produces polynomial of second order degree for Eq. (8), that can be resolved directly. It is important to note that, the stoichiometry of the system is always conserved. Moreover the solution of the equation never gives negative concentrations due to the limiting conditions (see Eqs. (9), (10)). Non physical solutions have never been met, but round off errors can give a null value of x. This problem has been avoided calculating approximate solution for x. For a second degree equation ax 2 + bx + c = 0, if a · c b2 , below the round off error, the direct solutions of the equation are x1 = 0 and x2 = −b/a, as in the case with c = 0. This problem can be easily solved remembering that x1 · x2 = c/a. Because x2 can be considered precise to the last digit, x1 can be approximated as x1 = c/a/x2 = −c/b which correspond to the solution obtained neglecting the terms with the power of x higher than one. This result can be easily extended to polynomial of degree higher than two. 2.4. Multi-reaction systems The calculation of equilibrium for multi-reaction systems could start from the equilibrium of a single reaction. The method here proposed consists in calculating the equilibrium concentrations solving hierarchically one-reaction equilibrium problem at a time, until the convergence is reached. Probably, the algorithm converges independently on the sequence of reactions chosen in the calculation [3], but it can be optimized choosing an appropriate sequence. In most cases, some reactions can be neglected, because they are almost ineffective for the equilibrium composition. At this purpose, the reaction distance from the equilibrium can be a suitable parameter to choose the reaction order: higher the distance, further the reaction from the equilibrium. The algorithm for the solution of the equilibrium problem for multireaction systems can be organized in the following way: starting from an initial value of concentration, we calculate the reaction position (ϑ) and distance (δ) for all the reactions. Than
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
we calculate the equilibrium composition for the reaction with the maximum distance. We iterate the previous procedure starting from the last calculated concentrations values. The solution of the single reaction equilibrium must be calculated with a good accuracy, otherwise, due to the large differences in the equilibrium constants, the reaction just considered could be solved again uselessly, reducing efficiency and accuracy of the overall equilibrium calculation. As stated before, there is more than one possible reaction position. We have chosen two different definitions. They behave differently and present advantages and drawbacks. The first definition (ϑ1 ) is ϑ1j = Fj (0) and the second (ϑ2 ) Kj · Πrj − 1 Kj · Πrj < Πpj , Πpj Πpj ϑ2j = Kj · Πrj > Πpj , 1− Kj · Πrj 0 Kj · Πrj = Πpj .
(11)
(12)
Both the definitions have the advantage to be continuous functions of the extent of reaction and finite real numbers even if some reactants or products have zero concentration. Anyway they behave differently: ϑ1 is not limited, while ϑ2 is between −1 (when at least one reactant has vanishing concentration) and 1 (when at least one product has vanishing concentration). For this reason the two positions are useful in different situations. The following example will clarify the problem: consider a diatomic molecule A2 that can be dissociated (a: A2 → 2A) and ionized (b: − A2 → A+ 2 + e ). Starting with the condition [A2 ]0 = 1 + and [A2 ]0 = [e− ] = [A] = 0, we have that the first reaction positions are ϑ1a = Ka and ϑ1b = Kb and the second are ϑ2a = ϑ2b = 1. Therefore ϑ2 cannot be used to discriminate the two reactions and therefore the choice can be made only on random bases. Suppose, for example, that the A+ 2 concentration is very small compared to the atom concentration and that the first reaction chosen is the ionization process. What happens is that after the ionization, the dissociation reaction looses some molecules, changing the ionization equilibrium toward the recombination weakly, so that the convergence is very slow. On the contrary ϑ1 chooses the reaction with higher equilibrium constant, which is the reaction furthest from the equilib-
181
rium. The drawbacks of using ϑ1 are: the round off errors of reactions with high equilibrium constants can be higher than the value of ϑ1 for reactions very far from equilibrium with small equilibrium constants and therefore such reactions will never be considered. As a consequence, the solution obtained gives the right concentration only for majority species, while secondary species present large relative errors. To take advantage of both distances, we build up two routines, with the same algorithm but using different distances. Firstly we call the routine that uses ϑ1 till the convergence is reached. At this step, reactions with small equilibrium constant have not been considered. Therefore we call the routine with the position ϑ2 . This distance does not discriminate the reactions by the equilibrium constant, therefore the reactions neglected in the first step, because of the small equilibrium constant, are now calculated determining the species with low concentrations. Obviously, all the reactions are included in the calculation for both distances; so the algorithm decides the reactions to consider. The calculation for each distance should terminate when all the reaction distances are below the tolerance required. While for ϑ2 it is trivial to satisfy this condition, for ϑ1 some troubles arise, because it is difficult to define an adequate tolerance value for all the reactions due to the large differences in the equilibrium constants. At this purpose an additional convergence criterion can be adopted based on the relative concentration error. Moreover, it is possible to decide if ϑ1 has reached the precision limit by considering ϑ2 which is a relative evaluation. This means that the in the first step, when the reaction position is ϑ1 , we calculate also ϑ2 ; if it is smaller than a given tolerance, ϑ1 is of the order of the round off error and then must be considered null. 2.5. Addition of new reactions Even if the double distance approach speeds up the calculation of equilibrium composition, some problems cannot be avoided: for mixtures with many species and many reactions, some numerical traps can be found and as a consequence the convergence slows down. This is the case in which two species concentration (A, B in Eq. (13)) are coupled by successive reactions passing by an intermediate species (X) with
182
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
much lower concentration Ra : Rb :
A + ··· ⇒ X + ···, X + ··· ⇒ B + ···.
(13)
In these conditions each reaction produces very small changes with respect to the absolute variations of A and B species, and therefore the two reactions in Eq. (13) are repeated many times one after the other till the global transfer is obtained. As an example, the dissociation/ionization described in the previous section can be considered, where A stands for the atom, B for the molecular ion and X for the molecule. The solution of the problems can be obtained considering a new reaction as the combination of the two written in Eq. (13) obtained eliminating the common element X as in Ra,b :
A + ··· ⇒ B + ···.
(14)
Villars [3] already pointed out this problem, suggesting as solution the addition of a new reaction. This criterion has not been generalized and is based on the experience of the researcher. For this reason the Villars approach has been abandoned. In this paper we propose an automatic method to add as many reactions as needed with the result to improve considerably the convergence. It must be noted that, differently by matrix methods, hierarchical approach can use any number of reactions. To choose the reactions to combine we count the number of time each reaction has been solved after a given number of iteration. If one of these reactions is higher than a threshold, it is combined with the reaction with the higher number of repetition which matches the requirements in Eq. (13). In practice, adding a new reaction is finding a shortcut to the path that the system must follow to reach the equilibrium. The complete algorithm for calculating the equilibrium with the addition of new reactions works in this way (see the flow chart in Fig. 1). We have two nested loops. The external loop (a “while” loop) checks only if the convergence is reached, the internal loop (a “for” loop) repeats the calculation of single reaction equilibrium a given number of iterations (Nint ). Within the internal loop we count the number of time each reaction has been chosen. If a given number of repetition (Nrep ) is reached by a reaction, we check between the other reactions with the higher repetition number if they fulfill the requirement in Eq. (13) so that they
can be combined as described in Eq. (14). If two reactions are combined in a new one, we solve this new reaction and we reset the reaction counters of the combined reactions. When the internal loop is completed, we check if the convergence requirements are reached and we reset the counters for all reactions. The values chosen for Nint and Nrep weakly affect the convergence behaviors, but they must clearly fulfill the requirement Nint /Nrep 2, otherwise the condition to build the new reaction will never be met. As stated before, we decide to build the new reaction when two reactions are repeated one after the other as described in the dissociation–ionization example. This is not the only situation can give problems. A more longer series of reactions can be looped together giving similar troubles. The probability of finding this situation decreases while the number of the reaction series increases. The kind of series we can get with the nested loop technique depend on the ratio Nint /Nrep . 2.6. Constraints A short comment about constraints is needed. Constant pressure and constant volume constraints, or any other referring to the number of molecules in the system, can be easily included in the code introducing an appropriate normalization condition after each equation has reached the equilibrium. For example, for an ideal gas at constant pressure P , the partial pressures are normalized by the factor P / i Pi , while the constant volume condition does not require any normalization. For non ideal gases such as plasma, the equation of state must be corrected according to the Debye–Huckel (DH) theory [11] 1 = Pid + PDH , P = Kb T n − (15) 24πλ3d where λd is the Debye length, n the total number density of particles and Kb the Boltzmann constant. In the frame of this theory, the correction to the pressure depends on the charged particle density. As for the number of molecule constraint, also for the DH plasma, it is sufficient to write an appropriate normalization condition. For example, for the constant pressure DH plasma we should normalize the partial pressures by a
factor (P − PDH )/ i Pi . The convergence of the DH corrections is faster than equilibrium, because it is important only if the ionization degree is larger than or
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
Fig. 1. Flow chart of the equilibrium calculation algorithm.
183
184
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
close to one and this algorithm determine the majority species concentration very fast and accurately. We must say that the DH theory modify the equilibrium constants as well. This is not a problem in the model because the convergence of the DH correction is very fast and not very sensitive to charge concentrations. We must point out that the DH equation of state is in the perturbation theory framework, which means that the deviation from the ideal gas is small. We do not apply the algorithm to equations of state with strong perturbations. The constant temperature constraint obviously gives no problem while the constant energy one introduces a further complication due to the dependence of the equilibrium constants on the temperature. This problem can be treated in different ways but we do not have a definitive answer on it. The temperature can be updated after all the equilibrium has been calculated or during the calculation of equilibrium after each reaction or a given number of reactions has been solved. We choose this second approach, but we did not test the convergence properties in one case or the other.
putting in evidence the improvement due to the added reactions. At this purpose we have applied the method to ionized gaseous mixtures in the temperature range 1000–30000 K. In particular we have considered a pure nitrogen mixture (see the first three reactions in Table 1 and Fig. 2) and the air mixture (see Table 1 and Fig. 3). The first system is very simple, therefore we can investigate the behavior of the algorithm step by step, being each step (iteration) the resolution of a single chemical reaction. It must be observed that the results obtained for this system do not include multiple ionized particles even if at high temperatures they
3. Results and discussion In this section we will discuss the results obtained by the method described above in different cases,
Fig. 2. Equilibrium composition of a nitrogen mixture as a function of temperature.
Table 1 Equilibrium constants of the reactions considered at different temperature. Species considered: N2 , N, N+2 , N+ , N++ , N+++ , N++++ , O2 , O+2 , − + ++ +++ , O++++ , NO, NO+ , e− O− 2 , O, O , O , O , O Reactions
log(Keq )
(16)
103 K
5 · 103 K
104 K
2 · 104 K
3 · 104 K
N2 2N N2 N2 + e− N N+ + e− N+ N++ + e− N++ N+++ + e− N+++ N++++ + e− O2 O+2 + e− O2 + e− O− 2 O2 2O + − OO +e O+ O++ + e− O++ O+++ + e− O+++ O++++ + e− N2 + O2 2NO NO NO+ + e− O + e− O−
−85.93590 −165.3994 −153.5874 −329.5341 −537.5887 −883.7819 −125.6238 −7.681201 −32.49372 −144.1756 −392.8815 −623.5494 −885.2708 −18.65947 −94.88367 3.713479
5.036874 −16.59715 −14.39963 −50.56606 −93.27721 −160.5291 −9.595570 −16.56471 15.62565 −13.81455 −62.25641 −109.3668 −162.8089 −1.346583 −4.713501 −15.26020
16.94815 3.624266 4.127968 −14.51822 −36.50893 −68.89310 5.960819 −19.04409 21.72691 3.714115 −19.72256 −43.89081 −71.27454 0.6411276 7.757746 −18.91194
24.26754 14.49052 13.53105 4.296274 −7.222843 −22.26790 14.84219 −21.37386 25.42207 13.35034 2.214520 −10.36335 −24.62963 1.521508 15.15033 −21.88424
29.62973 18.13983 15.70289 10.86431 3.104701 −6.460977 18.37340 −22.68414 28.97627 16.15614 9.825148 1.235014 −8.567586 1.581945 18.43401 −24.35636
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
185
Fig. 3. Equilibrium composition of air mixture as a function of temperature.
become important. In the second one we want to investigate the behavior of the algorithm for realistic large systems. In Table 1 we have reported also the equilibrium constants for each reaction at different temperatures showing the stiffness of the equilibrium problem. The calculation of the equilibrium constants is made starting from energy levels and partition functions. There are different aspects that need to be verified. In particular we consider the double distance approach and the addition of a new reaction. We consider four cases which differ by the reaction position considered, discussing the convergence: (1) the distance is ϑ2 . The first reaction solved is molecular ionization (reaction 2 in Table 1); (2) the distance is ϑ2 . The first reaction solved is the molecular dissociation (reaction 1 in Table 1); (3) the distance is ϑ1 . (4) combined solution of ϑ1 and ϑ2 .
The evolution of the molar fractions toward the convergence values with the number of reactions solved for pure N2 mixture has been reported in Figs. 4–8 at gas temperature T = 30000 K. As initial composition molecular nitrogen at p = 1 atm has been considered, which is very far from the equilibrium composition. The equilibrium composition has been calculated at constant volume. In Figs. 4a, b it is reported the gas composition as a function of the number of reactions solved in case 1. To reach the convergence (gray lines) all the species needs more than 10000 iterations. The N2 concentration (Fig. 4b) presents strong oscillations. This is exactly equivalent to the example of dissociation/ionization system discussed in Section 2.4. In this case the first reaction solved is the N2 ionization. As a result, the gas is made mainly by ionized molecules. Therefore the second reaction solved is the dissociation. These steps push the solution farther from the equilibrium, being N+ the dominant species. For
186
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
(a)
(b)
Fig. 4. Evolution of the molar fractions with the number of reactions solved for N2 system at T = 30000 K, case 1. The gray lines represent the equilibrium values at convergence.
(a)
(b)
Fig. 5. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 2. The gray lines represent the equilibrium values at convergence.
this reason the system goes back to the equilibrium changing slowly the ionized molecule and atomic nitrogen concentrations, alternating the solution of N2 ionization and dissociation. Each couple of reactions changes slowly the concentrations (N and N+2 ) be-
cause the intermediate species (N2 ) is negligible. After 10000 reactions solved the system is still very far from equilibrium and ionization and dissociation reactions are still solved alternatively making N2 concentration to oscillate on an interval of some orders of magnitude.
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
(a)
187
(b)
Fig. 6. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 3. The gray lines represent the equilibrium values at convergence.
The convergence changes if the first reaction solved is the dissociation (case 2), as can be observed in Figs. 5a, b. The majority species (N+ , N) reach the equilibrium values after few iterations (3 or 4), while minority species (N2 , N+2 ) need more than 1000 steps to converge. Also in this case there are strong oscillations of the N2 convergence series due to the alternate solution of molecule dissociation and ionization as in case 1. These results show the limit of ϑ2 : it is not appropriate to decide which is the most important reaction. If the right choice is made at the first step, the convergence can be speeded up of one order of magnitude. This problem disappears using ϑ1 (case 3, see Figs. 6a, b). In fact, at the beginning, the reaction farthest from the equilibrium is the one with highest equilibrium constant, in this case the dissociation. By using ϑ1 , the convergence values of majority species are reached after few iterations while minority species concentrations converge slowly (more than 1000 iterations), showing the same behaviors as case 2. This is a very important result, because, using ϑ1 , if only high concentration species are relevant, the convergence can be reached rapidly. In the last case (4, see Fig. 7) using ϑ1 the system needs few iterations to reach the convergence, but mi-
Fig. 7. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 4. The gray lines represent the equilibrium values at convergence.
nority species present wrong concentrations. On the other hand, using ϑ2 , also minority species concentrations are calculated with the required tolerance, even if they are very small. The oscillations in N2 concentration observed in all the cases discussed above put in evidence the presence
188
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
of some troubles in this algorithm. Adding new reactions introduces a shortcut to the solution, eliminating the oscillating behavior. The set of reactions used in the previous section is minimal, i.e. any further reaction is a linear combination of the minimal set. This approach is necessary for standard methods which needs to solve linear systems while non-minimal set leads to overdimensioned systems. In our algorithm only one equation at a time is solved, therefore there is no problem in using non-minimal sets. The problem is then to choose the set of reactions which leads to the equilibrium solution as fast as possible but this
Fig. 8. Evolution of the molar fractions with the number of reactions solved for the pure N2 system at T = 30000 K, case 1 with the addition of new reactions. The gray lines represent the equilibrium values at convergence.
choice is not possible a priori, because it depends on the thermodynamic conditions. In Fig. 8 we have reported the results obtained in case 1, adding new reactions automatically. In about 20 reactions solved, the equilibrium compositions are reached for all the species. It must be noted that there are only few oscillations in N2 concentrations, differently from Fig. 4b, because with the aid of the added reactions the concentrations jump directly to the correct values. In Tables 2, 3 we have reported the number of reactions to reach the convergence for a required tolerance (10−13 ) in the four cases listed above with (+R) and without new reactions for different temperatures,
Fig. 9. Evolution of the absolute value of ϑ2 with the number of iterations for air mixture at 1 atm and T = 10000 K, in different cases.
Table 2 Number of reactions solved to reach the convergence in different cases (see Section 3) and for different temperatures with (+R) and without the addition of further reactions for the pure nitrogen mixture. The tolerance required is 10−13 Cases
T = 103 K
T = 5 · 103 K
T = 104 K
T = 2 · 104 K
T = 3 · 104 K
1 1+R
3 3
31 36
31 27
1396 35
46506 24
2 2+R
3 3
30 28
29 33
818 31
24066 19
3 3+R
1 1
37 41
27 35
564 15
9328 12
4 4+R
3 3
31 57
27 66
865 52
25851 22
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
189
Table 3 Number of reactions solved to reach the convergence in different cases (see Section 3) and for different temperatures with (+R) and without the addition of further reactions for the air mixture. The tolerance required is 10−13 Cases
T = 103 K
T = 5 · 103 K
T = 104 K
T = 2 · 104 K
T = 3 · 104 K
1 1+R
59 59
510 302
2768 294
22591 684
>200000 435
2 2+R
57 62
520 406
2592 309
14589 833
>200000 327
3 3+R
3 3
641 128
2832 414
1277 171
19038 107
4 4+R
51 51
467 172
2394 328
13926 393
>200000 337
where the plasma is in completely different conditions, respectively for the pure N2 and air mixtures. In Fig. 9 we have reported the evolution of the absolute value of ϑ2 with the number of iterations for air mixture at 1 atm and T = 10000 K, in the cases labeled 1, 4 and 4 + R. The decreasing of |ϑ2 | is exponential and moreover the addiction of new reactions speeds up the convergence rate. In this particular case, the case 1 + R is faster than 4 + R at low tolerance. This is not a general case as can be observed in Tables 2, 3. The main advantage of the automatic addition of new reactions can be observed when the equilibrium composition is far from the initial values, reducing the number of iterations of many orders of magnitude. It must be noted that independently on the case considered, the addition of new reactions retains the number of iteration in the same order of magnitude. It is worth noticing that using ϑ1 the algorithm converges faster. The drawback is that very low concentrations are miscalculated. Nevertheless in some applications, such as equilibrium fluid dynamics or transport properties calculations, this approach can be used because only majority species are important.
4. Conclusions In this paper we have presented an efficient stoichiometric, hierarchical method to calculate single phase equilibrium composition. The novelty of this method respect to previous hierarchical algorithms is the automatic addition of new reactions in order to improve the speed of convergence. We have studied two
gas mixtures in different conditions to validate the algorithm, analyzing the differences with traditional approaches. The results obtained show that the algorithm is very efficient in all the conditions, independently on the initial values. We have considered ideal gas approximation but the method has been also applied for non ideal plasma with the Debye–Huckel correction. The method is very suitable for non ideal plasma calculation where the correction terms depend on the free electron concentration. The automatic addition of new reactions makes the code very efficient in a wide range of conditions. We have used it to calculate air properties [12] in the temperature range 50–105 K and in the pressure range 10−3 –103 atm. In all these cases the code never fails differently from STANJAN [13] and NASA [14] whose convergence is related to the constraints used [2]. Moreover the convergence of the method is enforced by the fact that it never calculates at any iterations non physical solutions such as negative concentrations. This code is the first attempt after about 50 years [3] to propose a modern version of a hierarchical approach that is automatic and fast enough to be of practical use. Recent efforts on the chemical equilibrium calculation are devoted to multiphase systems [7–10]. All these methods are based on the simultaneously solution of the non linear equation mainly considering Gibbs energy minimization. In our opinion the method here presented can be easily extended to multiphase systems.
190
G. Colonna, A. D’Angola / Computer Physics Communications 163 (2004) 177–190
References [1] W.R. Smith, R.W. Missen, Chemical Reaction Equilibrium Analysis: Theory and Algorithms, John Wiley and Sons, New York, 1982. [2] P.S. Bishnuy, D. Hamirouney, M. Metghalchiy, J.C. Keckz, Constrained-equilibrium calculations for chemical systems subject to generalized linear constraints using the NASA and STANJAN equilibrium programs, Combust. Theory Modelling 1 (1997) 295. [3] D.S. Villars, in: G. Bahn, E. Zuckowsky (Eds.), Proceedings of the First Conference on the Kinetics Equilibria and Performance of High Temperature Systems, Western States Section of the Combustion Institute, Butterworths Scientific Publications, Washington DC, 1960. [4] D.S. Villars, A method of successive approximations for computing combustion equilibria on high speed digital computer, J. Phys. Chem. 60 (1959) 521. [5] D.R. Cruise, Notes on rapid computation of chemical equilibria, J. Phys. Chem. 68 (1964) 3797. [6] K. Meintjes, A.P. Morgan, Performance of algorithms for calculating the equilibrium composition of a mixture of gases, J. Comp. Phys. 60 (1985) 219.
[7] C.M. McDonald, C.A. Floudas, GLOPEQ: A new computational tool for the phase and chemical equilibrium problem, Comput. Chem. Engrg. 21 (1997) 1. [8] Y. Sofyan, A.J. Ghajar, K.A.M. Gasem, Multiphase equilibrium calculations using Gibbs minimization techniques, Ind. Eng. Chem. Res. 42 (2003) 3786. [9] A.V. Phoenix, R.A. Heidemann, A non-ideal multiphase chemical equilibrium algorithm, Fluid Phase Equil. 150–151 (1998) 255. [10] R.N. Mioshi, C.L. do Lago, An equilibrium simulator for multiphase equilibria based on the extent of reaction and Newton– Raphson method with globally convergent strategy (SEQEx2), Anal. Chim. Acta 334 (1996) 271. [11] L.D. Landau, Statistical Physics, Elsevier Science, 1979, Chapter 7. [12] M. Capitelli, G. Colonna, C. Gorse, A. D’Angola, Transport properties of high temperature air in local thermodynamic equilibrium, Eur. Phys. J. D 11 (2000) 279. [13] W.C. Reynolds, The element potential method for chemical equilibrium analysis: implementation in the interactive program STANJAN, Stanford University Report ME 270 HO no 7, 1986. [14] S. Gordon, B.J. McBride, Computer program for calculation of complex chemical equilibrium compositions and applications, NASA Reference Publication 1311 (1994).