DYNAMIC APPROACH TO THERMODYNAMIC SIMULATION OF CHEMICAL SYSTEMS:

FROM TRUE EQUILIBRIUM TO TRUE CHAOS.

B. Zilbergleyt,

System Dynamics Research Foundation, Chicago

E-mail:

ABSTRACT.The paper presents new model of equilibrium in open chemical systems suggesting linear dependence of the reaction shift from true equilibrium on the external thermodynamic force. Basic equation of the model includes traditional logarithmic term and non-traditional parabolic term. Behavior of open system is defined by relative contributions of both terms. If classical logarithmic term prevails open system can be in equilibrium. Increased weight of non-classical chaotic term leads open system to bifurcations and chaos. Area of open equilibrium serves as a water shed between true equilibrium and bifurcations and eventual chaos. In simple isolated equilibrium the chaotic term equals to zero turning the equation to traditional form of constant equation. Formally this term matches excessive thermodynamic function showing linear relationship between reaction extent in open equilibrium and logarithm of coefficient of thermodynamic activity. In form of chaotic term it reveals new behavioral features of open systems while method of coefficients of thermodynamic activities hides them to keep open system disguised as isolated entity. Discovered relationship prompts us to use the coefficient of linearity in combination with reaction shift rather than activity coefficients. This coefficient in open equilibrium can be calculated in a simple way using basic equation and reaction shift as independent variable. Numerical data obtained by various simulation techniques have proved premise of the new theory and following from it method of chemical dynamics.

INTRODUCTION: BACK TO CHEMICAL DYNAMICS.

Nowadays we know that chemical self-organization happens in a vaguely defined area “far-from-equilibrium”[1], while classical thermodynamics defines what is frozen at the point of true equilibrium. What occurs in between?

“True”, or “internal” thermodynamic equilibrium is defined by current thermodynamic paradigm only for isolated systems. That’s why applications to real systems often lead to severe misinterpretation of their status, bringing approximate rather than precise results. A few questions arise in this relation. Is it possible to expand the idea of thermodynamic equilibrium to open systems? How to describe and simulate open equilibrium in chemical systems? Is there any relationship between deviation of a chemical system from “true” equilibrium and parameters of its non-ideality?

Current paradigm of chemical thermodynamics handles open chemical systems employing coefficients of thermodynamics activities. The coefficients as well as the whole concept of equilibrium, based on Guldberg-Waage’s equation, may be derived from the probability theory [2]. Indeed, chemical system with only one type of collision is a simple chemical system: one type of collision – one reaction – one outcome, and chance of the reaction to happen is proportional to P(A) - probability of reactants to collide, equal to the particles mole fractions product. In complex systems with multiple interactions, according to the Bayes’ theorem, instead of P(A) we consider P(A/B) – conditional probabilityof collision A given another collision, B, took place. The ratio=P(A/B)/P(A) defines a coefficient, currently known as coefficient of thermodynamic activity.

Regardless of the initial idea, introduction of coefficients of thermodynamic activity was and still is a great contribution to chemical thermodynamics mainly because it allowed scientists to keep expressions for G and constant equation unchanged still holding open systems in a mask of isolated entities. The name of the Law of Active Masses became more clear -  may be considered a fraction of a substance available to participate in chemical reaction. Later on value of  was tied to excessive thermodynamic functions though experiment is often the best way to obtain it.

The values to be introduced and discussed in this paper are sensitive to isolation or openness of the systems where the chemical processes are to run. An isolated system will be referred to as system with possible true thermodynamic equilibrium. Simple chemical system allows only one chemical reaction to run towards true thermodynamic equilibrium. Thermodynamic state of such a system is defined by two thermodynamic (e.g., P, T) and one chemical (reaction coordinate parameters. For the sake of clearness, simplicity and to keep it closer to reality, open system will be considered also allowing for one chemical reaction to run thus being a simplesubsystem – a part of bigger and more complex chemical system. The system contains a set of similar entities as subsystems and constituted by chemical (or thermodynamic) interaction between them. For example, in some cases it is very convenient if reaction of formation of one chemical substance from elements runs in each subsystem. This way number of subsystems will equal to the number of chemical species in the system. We define thermodynamic equilibrium in open system as open equilibrium. It is the central idea of the present paper, and we will get its basic equations, show where it leads the open system and how to bring the ideas to numbers.

To do so we will use currently almost neglected method de Donder that introduced thermodynamic affinity interpreting it as a thermodynamic force and considering the reaction extent a “chemical distance” [2]. It is more convenient to use redefined reaction coordinate dj=dnkj/kj, instead of dj=dnkj/kj by de Donder, or reaction extentj =nkj/kj in increments. Value of nkj equals to amount of moles, consumed or appeared in j-reaction between its two arbitrary states, one of them may be the initial state. The ij value equals to a number of moles of k-component, consumed or appeared in an isolated j-reaction on its way from initial state to true equilibrium and may be considered a thermodynamic equivalent of chemical transformation. This value is unambiguously related to standard change of free Gibbs’ energy or equilibrium constant of reaction and represents natural and the only result of thermodynamic simulation of the reaction equilibrium. Above redefined value of the reaction extent remains the same being calculated for any component of a simple chemical reaction; the only (and easily achievable by appropriate choice of the basis of the chemical system) condition for this is that each chemical element is involved in only one substance on each side of the reaction equation like in reactions of formation from elements. Now, in our definition j is a dimensionless chemical distance (“cd”) between initial and running states of j-reaction, 0≤j≤1, and thermodynamic affinity A= (G/)p,Tturns into a classical force by definition, customary in physics and related sciences.

Chemical reaction in isolated system is driven only by internal force (eugenaffinity, Aij). True thermodynamic equilibrium occurs at A´ij = 0 and at this point ´j = 1. Reactions in open system are driven by both internal and external (Aej ) forces. The external force originates from chemical or, in general, thermodynamic (also due to heat exchange, pressure, etc., affecting thermodynamic parameters of its state) interaction of the open system with its environment. Linear constitutional equations of non-equilibrium thermodynamics at zero flow give us a condition of open equilibrium with resultant affinity

A*ij + ie A*ej = 0, (1)

where ie = aii /aie is a ratio of the Onsager/Kasimir coefficients [3]. The accent mark and asterisk relate values to isolated (“true”) and open equilibrium correspondingly, and indices ‘i’ and ‘e’ define internal and external variables and functions correspondingly.

In this work we will use only one assumption which in fact slightly extends the hypothesis of linearity. Given a relation between the reaction shift from equilibrium j =1j and external thermodynamic force causing this shift, we will supposeat the first approximation that the reaction shift in the vicinity of true thermodynamic equilibrium is linearly related to the shifting force

j = ie Aej . (2)

Recalling that Ai = (Gi /i ), or Ai = (Gi /i ) and substituting (2) into (1), we will get an intermediate expression

– (G*ij /*j) + (ie/ie)*j = 0. (3)

After a simple transformation, retaining in writing only j for j and j for j we turn it into

G*ij + bie *j *j = 0, (4)

where bie = ie /ie). As usually, G*ij = G0ij + RTln*j (, *j) , and corresponding constant equation is

 RTlnKj + RTln*j (, *j) + bie *j * j = 0, (5)

where *j (, *j) is the product of mole fractions in open equilibrium expressed via j-reaction extent. Please notice that equilibrium constant Kj = 'j (, ).

Obviously product *j *j is dimensionless, therefore bie has dimension of energy. To bring (5) to more symmetric form we introduce a new value - the “alternative” temperature of the open system

Ta = bie / R, (6)

R is universal gas constant. The value of Ta is introduced in this work for convenience and symmetry; we cannot clearly explain its physical meaning at this moment. The logarithmic term contains traditional thermodynamic temperature Tt , and (5) turns to

RTt ln`j + RTt ln*j + RTa *j*j= 0.

Now, dividing (7) by (RTt ), and denoting = Ta /Tt, we transform equation (7) into

ln [`kj,1)/ kj , *j)]  j*j*j = 0. (8)

______

 “Vicinity” in this case is certainly not less vague than “far-from-equilibrium”. Some discussion see below.

So, as soon as chemical system becomes open, appropriate constant equation includes a non-linear, non-classical term originated due to system’s interaction with its environment.

What opens up immediately is a similarity between the non-classical term of (5) and the well known product rx(1x) from so called logistic map [4], which is one of the most convincing equations leading to bifurcations and chaos. We refer alternative temperature to as chaotic temperatureTch and j as reduced chaotic temperature.

Being divided by *j , this equation expresses linearity between the thermodynamic force and reaction shift

{ln [`kj,1)/ kj , *j)]}/*j =  j*j, (9)

both parts of it are new expressions for the thermodynamic force. Equality (9) contains condition of open equilibrium as balance of internal and external thermodynamic forces. Equation (8) can be written in more general form as

*= *(10)

It is easy to see that in case of isolated system *= external thermodynamic force equals to zero, and (8) turns to the normal constant equation. For better understanding of internal relations between (9) and constant equation one should recall once again that serving as a parameter in(10), is the only output from the solution to constant equation (because n` n0, where right side contains equilibrium and initial mole amounts).

INVESTIGATION OF THE FORCE-SHIFT RELATIONSHIP.

First, consider the force expression from equation (9). Its numerator is a logarithm of a combination of molar

fraction products for a given stoichiometric equation. The expression under the logarithm is the molar fraction product for ideal system divided by the same product where kj replaced by *jkj due to the system’s shift from “true” equilibrium. Table 1 represents expressions for thermodynamic forces of some simple chemical reactions with initial amounts of reactants A and B both equal to one mole. Graphs of the reaction shifts vs thermodynamic forces are shown at Fig.1. One can see visually distinctive linearity (actually, quasi-linearity) on shift-force curves. Extent of the linearity region depends on the value.

Going down to real objects, consider a model system containing a double oxide MeORO and an independent reactant I (for instance, sulfur) such that I reacts only with MeO, while RO restricts reaction ability of MeO

Table 1.Thermodynamic forces {ln[`kj,1)/ kj , *j)]}/ (eq. 8) for some simple chemical reactions.

Initial amounts of reactants are taken equal to 1 mole and initial products to zero for simplicity.

Reaction equation. / Thermodynamic force
A + B = C / {ln {(1/(2)/(2)][(1)/(1)]2}}/
A + 2B = C / {ln{(1/(1)/(1 [(12)/(12)]2}}/
2A + 2B = C / {ln{(1/[(23233[(124}}/

Fig. 1. Shift of some simple chemical reactions from true equilibrium  vs. shifting force F, kJ/m·cd. Reactions, left to right, values of in brackets(curves follow right to left) : A+B=AB (, 0.9), A+2B=AB2 ( 2A+2B=A2B2 (. Also, linear areas on the curves give an estimation of how far the “vicinity of equilibrium” extents.

and frees in the reaction as far as MeO is consumed. Two competing processes are in equilibrium in the isolated system - decomposition of MeORO, or restricting reaction: MeORO = MeO + RO, and leading reaction: MeO + I =*L, the right side in the last case represents a sum of products. Resulting reaction in the system obviously is MeORO + I =*L + R.

To obtain numbers for real species, we used thermodynamic simulation (HSC Chemistry for Windows) in the model set of substances. The Is were S, C, H2, and MeORO were double oxides with symbol Me standing for Co, Ni, Fe, Sr, Ca, Pb and Mn. As restricting parts RO were used oxides of Si, Ti, Cr, and some others. Chosen double oxides had relatively high negative standard change of Gibbs’ potential to provide negligible dissociation in absence of I. For more details see [5]. Some of the results for reactions (MeORO+S) are shown on Fig.2. In this group value of (G0C / *L)Iplays role of external thermodynamic force regarding the (MeO+S) reaction.

The most important is the fact that in both cases the data, showing the reality of linear relationship, have been

received using exclusively current formalism of chemical equilibrium where no such kind of relationship was ever assumed. It is quite obvious that linear dependence took place in some cases up to essential deviations from equilibrium. We call the method, described in the present work (also including original de Donder’s approach), a force-shift method for explicit usage of chemical forces, originally introduced as thermodynamic affinities, or a method of chemical dynamics (MCD). Results shown on Fig.1 and Fig.2 prove the premises and some conclusions of the theory.

Fig.2. Reaction shift * vs. force F= G0 MeO RO / *, kJ/m·cd, 298.15K, direct thermodynamic simulation.

Points on the graphs correspond to various RO. One can see a delay along x-axis for CaORO.

Within current paradigm of chemical thermodynamics, constant equation for non-ideal system with k1 is

G0j =  RTt ln *k RTtln*xk,

and xk are molar fractions, power values are omitted for simplicity. The non-linear term of the equation (8) also belongs to a non-ideal system, and comparison of (8) and (15) leads to following equality in open equilibrium with precision of the sign

j *j ln *k)/*j. (16)

If *k<1, the minus sign should be placed on the left side. This result is quite understandable. For instance, in case of MeORO the chemical bond between MeO and RO reduces reaction activity of MeO; the same result will be obtained for reaction (MeO + I) in absence of RO and with reduced coefficient of thermodynamic activity of MeO. Now, to avoid complexity and using only one common component MeO in both subsystems, the relationship between the shift of the (MeO +I) reaction and activity coefficient of MeO is very simple

* = (1/)(ln *)/*], (17)

where *, and* are related to the (MeO +I) reaction, and(ln *)/*] is external thermodynamic force acting against it (divided by RTt). This expression for the force as well as the total equation (17) are new. This equation connects values from the MCD with traditional values of chemical thermodynamics. Yet again, at *L= 0 we have immediately *=1, and vice versa, a correlation, providing an explicit and instant transition between open and isolated systems. In case of multiple interactions one should expect additivity of the shift increments, caused by interaction with different reaction subsystems, which follows the additivity of

appropriate logarithms of activity coefficients. This is also proven by simulation.

Data for Fig. 3 were obtained using two different methods of thermodynamic simulation. I-simulation relates to an isolated (MeORO+I) system with real MeO and RO and MeO= 1 in all cases. In O-simulation a ______

 The numerator is free energy of formation of double oxides MeORO from oxides MeO and RO.

combination of |MeO+Y2O3+I| represented the model of open system where RO was excluded and replaced by yttrium oxide, neutral to MeO and I to keep the same total amount of moles in the system as in I-simulation and avoid interaction between MeO and RO, which role was played by Y2O3. Binding of MeO into double

compounds with RO, resulting in reduced reaction ability of MeO was simulated varying . I-simulation provided a relationship in corresponding rows of the *L - * values, and O-simulation - with *L - G0MeO RO

Fig. 3. *vs. (-ln *) (I-simulation, x) and vs. (G0MeORO/*) (O-simulation, o), (MeOR+S). PbO and CoO at 298K, SrOat 798. Curve for SrO shows light delay along the x-axis.

correspondence. Standard change of Gibbs’ potentialG0MeO RO, determining strength of the MeORO bond, was considered an excessive thermodynamic function to the reaction (MeO+I).

We have calculated some numeric values from the data used for plotting Fig.2, they are shown in Table 2. It is worthy to mention that the range of activity coefficients, usable in equilibrium calculations, seems to be extendible down to unusually low powers (see Fig.1).

Strong relation between reaction shifts and activity coefficients means automatically strong relation between shifts and excessive thermodynamic functions, or external thermodynamic forces. Along with standard change of Gibbs’ potential we also tried two others – the Qwhich was calculated by equation (17) with *, used in the O- simulation, and another, G*MeORO, found as a difference between G0MeORO and equilibrium value of RTtln*xk. Referring to the same *, all three should be equal or close in values. Almost ideal match, illustrating this idea, was found in the CoORO - S system and is shown on Fig.5. In other systems all three were less but still are close enough. Analysis of the values, which may be used as possible excessive functions, shows that the open equilibrium may be defined using both external (like G0MeORO) and internal (the bound affinity, see [4]) values as well as, say, a neutral, or general value like a function calculated by (17) at given activity coefficient.

Table 2.

Reduced temperatures, standard deviations and coefficients of determination between *and (ln *) in

some MeORO-S systems. Initial reactants ratio S/MeO = 0.1.

CoO*RO / SrO*RO / PbO*RO
Tt , K / 298.15 / 798.15 / 298.15
 / 40.02 / 6.54 / 3.93
Standard deviation, % / 8.99 / 2.99 / 6.80
Coefficient of determination / 0.98 / 0.99 / 0.97

In principle all three may be used to calculate or evaluate L. This allows us to reword more explicitly the problem set in the beginning of this work and explain the alternative temperature more clear. It is easy to see that equation (17) represents another form of the shift-force linearity.

Recalling that  = (Tch /Tt), one can receive

* = [1/(RTch)]• (QE / *), (18)

where QE is a general symbol for excessive thermodynamic function. It means that the shifting force is unambiguously related to the excessive function, and the alternative temperature is just inverse to the coefficient of proportionality between the force and the shift it causes. The product RTch has dimension of energy while and  are dimensionless.

To compare values of received with different methodswe put them together in one Table 3. Abbreviations in the table means: R-simulation – thermodynamic simulation of homologous series of reactions with MeORO, (ElRea)-simulation – abstract simulation of elemental reactions with corresponding reaction equation and varied , and simulation stands for thermodynamic simulation with artificially varied coefficients of thermodynamic activity. In some cases (like in reaction 2CoO+S=2Co+SO2 at 298K) results match very well, but in some cases, like in the example below, the match is not very good.