Ind. Eng. Chem. Res. 2002, 41, 2931-2940
2931
Calculation of Entropy Production in Process Models Pertti S. Koukkari*,† and Simo S. Liukkonen‡ VTT Chemical Technology, P.O. Box 1404, Espoo FIN 02044-VTT, Finland, and Laboratory of Physical Chemistry and Electrochemistry, Helsinki University of Technology, P.O. Box 6100 FIN-02150 HUT, Finland
A thermochemical method is presented by which multiphase processes can be simulated with concurrent calculation of the Gibbs energy of the reactive mixture during a chemical change. Algorithmic constraints are set for the overall reaction kinetics when the Lagrange method of undetermined multipliers is used to minimize the Gibbs energy of the multicomponent system. Consequently, the chemical change is calculated as a series of successive “virtual” states, which follow the extent of the overall reaction. From the Gibbs energy of these intermediate states, other thermodynamic quantities for the changing chemical system can be deduced, and the entropy production of the process can be calculated. A particular process model was developed for two isothermal aqueous systems and for a nonisothermal high-temperature process. The entropy production in the reactive systems is presented, and the validity of the simulation models is assessed in terms of their time-dependent Gibbs energy and entropy profiles. Introduction The simulation of chemical processes and reactive flows by efficient computer programs has become everyday practice in both research and industry. For processes that involve time-dependent changes in both temperature and chemical composition, most often, models based on global reaction rates are used. The mechanistic approach, which involves elementary reaction kinetics, is also frequently applied. In these models, thermodynamic quantities are generally used to support material and energy balances and to check final equilibrium conversions. Obviously, no model should violate the conservation of energy and mass, and the chemical changes described in the models should proceed toward a defined state of equilibrium composition. The latter condition is satisfied provided that the time span of the chemical reactions is sufficient and that long-term metastable states are not present in the system under consideration. In this work, the thermochemical simulation method, previously used to combine dominant reaction rates with multicomponent thermodynamic Gibbs energy models,1 was applied to evaluate entropy production in three different processes. Two of the examined systems involved reactions in aqueous solutions, and the third was a high-temperature flame process. For each of these processes, a dynamic or stationary process model was produced by using Gibbs energy simulation combined with measured global reaction rates. The reaction rate constraints were applied to reactant species to enable the calculation of the intensive properties of the reaction mixture during the reaction path. The corresponding changes in Gibbs energy and entropy production were also calculated, and the behavior of these fundamental quantities during the multicomponent chemical change was evaluated against general thermodynamic principles. * To whom correspondence should be addressed. E-mail:
[email protected]. Tel.: +358-9-456 6366. Fax: +358-9456 6390. † VTT Chemical Technology. ‡ Helsinki University of Technology.
Thermochemical Simulation Method with the Extended Reactant Matrix. The thermochemical method, in general, involves the minimization of the Gibbs energy of a given thermodynamic system in terms of the molar amounts of its constituents.2-4 The temperature and pressure of the system must be defined, and the sole additional constraint needed is the elemental balance of the chemicals involved. The minimization of the Gibbs energy can be performed when the chemical potentials of the substances are known as functions of their standard-state and temperature-dependent excess Gibbs energy data. The Gibbs energy of a multicomponent system is written in terms of the chemical potentials as
G)
∑R ∑k nRk µRk
(1)
where µRk ) µRk (T, p, nRk ) is the chemical potential of the species k in phase R and nRk is its molar amount. The Gibbs energy is an extensive state variable, and the chemical potential µRk is the partial molar Gibbs energy of constituent k. The chemical potential is commonly written in terms of the standard molar enthalpies and entropies and temperature-dependent molar heat capacities of the constituents. This expression as input data for the Gibbs energy calculations provides the user not only with the Gibbs energy but also with other state functions of the system, such as the enthalpy. Therefore, a general advantage of the nonstoichiometric Gibbsian method in practical process simulation is that the calculation simultaneously deals with the complete heat and chemical balance of the multiphase and multicomponent chemical change. The condition for the conservation of the amounts of elements is φ
N
∑R ∑k aRkj nRk ) bj
(j ) 1, 2, ..., l)
(2)
where the coefficients aRkj are the elements of the stoichiometric matrix of the multicomponent system, φ is the number of phases present, and N is the total
10.1021/ie010498t CCC: $22.00 © 2002 American Chemical Society Published on Web 05/08/2002
2932
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002
number of constituents of the system. In general, the aRkj coefficients represent the chemical elements, but for a number of applications, it is also useful to introduce a chemical substance, or ion, as a basic component of the minimization procedure. Thus, l is the number of elements or, alternatively, the number of components of the minimization system. The Lagrange method of undetermined multipliers, as described more closely in Appendix A, is often used to minimize the Gibbs energy with the elemental balances as additional constraints. When the Gibbs energy minimization is performed by imposing the condition in eq 2, a composition corresponding to the global minimum in terms of the molar amounts of the constituents is produced. When the Lagrange method is used for the minimization, an extension of the stoichiometric matrix can be used to calculate a reaction path in a multicomponent system.5 In this method, the kinetically constrained reactants are introduced to the system as inert “image” constituents, which gives the extended mass balance condition φ N/
∑R ∑k aRkj nRk ) b/j
(j ) 1, 2, ..., l*)
(3)
where the increase in the total number of components from N to N* (N* > N) and in the total number of elements from l to l* (l* > l) forms the basis of an extended matrix. As the total mass mtot must be conserved, the additional condition between the global equilibrium system and the extended system is then l/
l
∑j
bjMj )
∑j b/j Mj
(4)
where Mj represents the molecular mass of the component j. The conversion of an image constituent to an active reactant is performed algorithmically in the minimization procedure.5 A similar balance for a more general matrix of a Gibbsian calculation can also be derived.6 The extended matrix allows for the conservation of both the mass and the thermodynamic properties of the system during the kinetically constrained Gibbs energy minimization. The essential feature of the method is that the species that is kinetically conserved is only allowed to change its molar amount (concentration) according to a given reaction rate. It is thus regarded as inert in the Gibbsian thermodynamic procedure.5,7 The reaction rates (Rr) are introduced as usual in terms of their rate coefficients, concentrations, and temperature
Rr ) F(kr, [A], [B], ..., [M])
(5)
The rate constant (kr) is a function of temperature and follows the Arrhenius equation
ln kr ) ln Rr -
βr T
(6)
where Rr and βr are parameters defining the frequency factor and the activation energy, respectively, of reaction r. The expressions for global reaction rate laws vary considerably according to the particular application under consideration. When each intermediate state during the reaction is calculated by Gibbs energy minimization, one should reach a reaction path for
Figure 1. Volume element of a chemically reactive system.
which the thermodynamic properties become defined. Consequently, rational behavior of the Gibbs energy and entropy in terms of the extent of the reaction is obtained. In Figure 1, a volume element of the dynamic chemical system is shown. It is assumed that no external force fields (such as an electric field) are present and that no external work is done for the system. The volume element is defined by its temperature, pressure, and chemical composition. The element encounters an entering reactant flow (m ˘ R), and the respective output of products (m ˘ P) is leaving the element (diffusive side streams are omitted for simplicity). The element is characterized by a measurable temperature and pressure, as well as by a composition in terms of the molar amounts of substances nk. Further, the exchange of heat (Q˙ ) between the volume element and its surroundings is shown. As the temperature, pressure, and molar amounts of substances are known, the Gibbs energy G ) G(T, p, nk) is defined for the volume element. The minimization of the Gibbs energy at given T and p in terms of the molar amounts is the traditional method for calculating the global equilibrium composition of the system. During the past few years, an increasing number of methods have emerged that apply the Gibbsian method for time-dependent process calculations. These include diffusion phenomena coupled with local equilibrium chemistry,8 a combination of flowsheet techniques with the Gibbs energy calculation,9 and a combination of reaction kinetic constraints with the Gibbsian multiphase procedure.1,5 In such calculations, the chemical process becomes described in terms of the subsequent minima of the Gibbs energy, taking into account the heat and mass transfer between the volume element and its surroundings, as well as the additional constraints due to the measured overall reaction rate. The optional constraints can then partly replace the elemental balance and T-p conditions of the original minimization procedure and are usually introduced with algorithmic techniques. The general procedure of the Gibbsian simulation is presented in Figure 2. From the Gibbs energy, it is then unambiguous to deduce other thermodynamic state properties of interest. If, during the simulation procedure, the volume element is considered as a function of time, the respective thermochemical properties also can be calculated as time-dependent functions. With the recent formalism of Kondepudi and Prigogine,10 the entropy production of an arbitrary isotropic system is presented as
S)-
∂G ∂T
dS ) diS + deS
(7) (8)
where the differentials diS and deS are related to the entropy production and entropy flow, respectively, of the thermodynamic system. For open systems, in terms of
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002 2933
Figure 3. Computed redox potential compared with measurements. The calculated pH and ionic strength (mol kg-1) are also shown.
Figure 2. Gibbsian method in process simulation.
the internal energy (U) and volume (V) of the system, one can write
deS )
dU + p dV T
diS ) -
1
-
1 T
N
∑1 µk denk
(9)
N
∑µk dink g 0 T 1
dnk ) dink + denk
(10)
Fe2+(/)(aq) f Fe2+(aq) f Fe3+(aq) + e-
(11)
Here, Fe2+(/) describes the image species, which is stoichiometrically isolated from the rest of the ferrousferric constituents of the system. The image has equivalent thermodynamic and activity data with the actual thermodynamic species [Fe2+], but because of the stoichiometric constraint, it remains in the thermodynamic calculation at its given concentration. The extended stoichiometric matrix used in the redox model with the appropriate image species Fe2+(/)(aq) is presented in Appendix B. During the discrete steps of the calculation, the image species is gradually changed to Fe2+, which, in the system of excess oxygen and sulfuric acid, is instantly oxidized to Fe3+. Thus, the reaction rate can be controlled by the algorithmic change Fe2+(/) f Fe2+ and the overall mass balance follows from eq 4. The rate law for the change is given by a second-order rate equation using the parameter values as described in Appendix B (see also, e.g., Iwai et al.13). The subsequent solution reactions in the multicomponent model are assumed to follow thermodynamic principles. Accordingly, the multicomponent model can be used to calculate the ionic activities of the redox couple during the extent of the reaction, as well as the changes in the ionic strength and pH of the solution. The Pitzer equation was used for the specific ion-ion interaction parameters in the activity calculation.14 In Figure 3, a typical model result is shown for the redox potential on the Ag/AgCl scale. The predicted values are compared with the redox potential actually measured at 5 bar (oxygen partial pressure). The computed values of pH and total ionic strength are also shown. The time-dependent redox potential calculated from the model is in general agreement with the measured points. The measured and calculated potential values, however, start to deviate from each other when the Fe3+ concentration increases. Two obvious reasons contribute to this trend: because of the lack of ion-ion interaction data, the model does not include specific ionic interaction parameters for the ferric ion,
and, for closed systems at constant temperature and pressure
deS )
dQ T
diS 1 dG )g 0 (T, p constant) dt T dt
terms of the changing activities of the ferrous-ferric couple in the sulfate system. The extent of the redox reaction was constrained in the model by using the measured reaction kinetics of the oxidation process. The ferrous ion was introduced to the calculation as a stoichiometrically isolated image species that was gradually allowed to react according to the scheme
(12) (13)
Correspondingly, the increment of the molar amount of substance k has also been divided to dink and denk. The latter of these is the increment that describes the exchange of substance k between the volume and its surroundings, with dink being the increment due to chemical changes within the volume element itself. The entropy production should remain greater than zero for any process without external work. For closed systems at constant temperature and pressure, the calculation of entropy production is particularly simple, as the time derivative of the Gibbs energy directly provides the desired quantity. In the absence of diffusive side streams (denk ) 0) the system can be treated as closed during any step of a sequential simulation and the entropy production in each sequence is then easily obtained by using eqs 8 and 12. Calculation Examples Dynamic Model of the Extent of a Redox Reaction in Aqueous Solution. The ferrous-ferric oxidation in aqueous sulfate systems takes place at moderate reaction rates when performed under pressurized autoclave conditions and at elevated temperatures. The features of the redox process are collected in Appendix B. In this study, the image component technique was used to calculate the change of the redox potential in
(14)
2934
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002
Thermodynamic Model for CaCO3 Dissolution. The thermochemical method was applied to model an aqueous system containing undissolved calcite which was then subsequently treated with small additions of 1 M sulfuric acid. The dynamic model of calcite dissolution was then designed according to the experimental results obtained. The applied thermodynamic model sequence was as follows: (1) The addition of sulfuric acid or carbon dioxide lowered the pH of the solution and brought about dissolution of the CaCO3 particulates. (2) The reactions in solution increased the pH. (3) The reactions in solution, including the release of gaseous CO2, were assumed to be fast when compared to dissolution. With this sequence, a simple reaction rate expression could be combined with the multiphase Gibbsian model to follow the effect of subsequent acidifications of the solution. The reaction kinetics of CaCO3 dissolution has been assessed by Compton and Daly.16 Considering different rate models for the solution process, they suggest that the flux of Ca2+ ions from the solid phase to the solution can be described by the simple rate equation
jCa2+ (mol cm-2 s-1) ) k1[H+] + k2 Figure 4. Gibbs energy, entropy, and entropy production during the redox process.
and basically for the same reason, the formation of sulfate complexes induced by the Fe3+ ion were neglected in the model. Regarding the simplicity of the multicomponent model, the result is satisfactory, as the reproducibility of the electrical state generally is tedious to maintain in the transient redox system. (In fact, as a result of the measurement series, it was concluded that, although the measured potential could be followed in terms of the extent of the redox reaction during single experiments, the reproducibility of the measured redox potential values was not sufficient to allow for a comparison of reaction rates under varying process conditions.11) The simple redox process takes place at constant temperature and follows the second-order rate equation in terms of its key reactant, the Fe2+ ion. The subsequent reactions in the solution can be assumed to be fast. From the thermochemical model, which applies the Gibbs function for all of the time-dependent intermediate states, one can deduce the values of both the Gibbs energy and the entropy during the redox process. The results are given in Figure 4, showing that the model is consistent with the requirement of having the Gibbs energy be a monotonically descending function in the constant temperature and -pressure processes. As a direct consequence, the entropy production of the redox reaction is positive and approaches zero when the reaction proceeds toward equilibrium. Dissolution of Calcium Carbonate by Sulfuric Acid. A dynamic model for the dissolution of calcium carbonate in sulfuric acid is of interest in the contemporary process of paper making from recycled pulp. The calcium carbonate is brought into the process with the recycled fiber, whereas the pH control of the process is still most often achieved with sulfuric acid and caustic soda. The reaction kinetics of the dissolution of CaCO3 was studied experimentally as described earlier.15 A brief summary of the experiments is given in Appendix C.
(15)
If the effect of the particle size is neglected, the increment of the dissolved Ca2+ ions becomes
d[Ca2+] ) (k1[H+] + k2)a dt
(16)
where [Ca2+] and [H+] are the concentrations of the Ca2+ and H+ ions, respectively; t is the time; k1 and k2 are the respective reaction rate parameters; and a is an additional parameter derived from the average particle size. In this form, the dissolution rate can be directly combined with the Gibbsian equilibrium model, which can be solved, e.g., by using discrete increments of soluble calcium carbonate. However, the image component method allows one to determine the contribution of the undissolved CaCO3 to the Gibbs energy of the system and, hence, makes possible the calculation of the thermodynamic functions of the intermediate states from the model. For the reactions in solution, no kinetic constraints are applied, and thus, the assumption in the dynamic multicomponent model is that the reactions in solution are at equilibrium and that the minimum of the Gibbs energy for the solution phase is attained. With this rather simple model, the pH variation caused by the subsequent acid additions can be simulated with adequate accuracy between the measured and calculated pH and Ca2+ concentration values (Figure 5). The use of thermodynamic states in the simulation is particularly useful to ensure the robust representation of the amplitude of the pH oscillations, ranging from pH 2.5 to pH 8 in the model. The increase of acid during different intervals brings the mixture to a low pH, which is determined by the amount of solution and the acid excess. Assuming fast acid dissociation as compared with calcite dissolution, the low pH is that of the dissolving solution, which also is the lowest attainable pH value. Correspondingly, in the upper end of the pH pulse, the simulated pH value slowly approaches the equilibrium pH, corresponding to the thermodynamic state where
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002 2935
Figure 5. Calculated pH and Ca2+ dissolution curves compared with measurements.
Figure 6. Calculated entropy production compared with pH curves of the calcite dissolution process.
the dissolving carbonate has entirely consumed the excess acid. The dots in the upper curve of Figure 5 indicate the concentration of Ca2+ ions in the solution as measured with an ion-selective electrode (ISE) and with atomic absorption spectroscopy (AAS) of microfiltered samples. The AAS results are more likely to represent the correct values, as the operation of the ISE was disturbed by the subsequent acidifications.15 The model is also thermodynamically consistent, as can be seen in Figure 6, as the entropy production in the constant-temperature system remains positive during the reaction periods and approaches zero when true equilibrium is attained. There is, however, some distortion in the subsequent entropy curves. A more specific kinetic model most likely would end up with more identical entropy production in the repeated acidification phases. Model of the Titanium(IV) Chloride Burner. Pure titanium dioxide (TiO2) is a bulk commodity that is used as a white pigment for paints, plastics, paper, and rubber. The particular reactor modeled in this work
was the titanium(IV) chloride oxidation burner. In the chloride process of TiO2 manufacture, rutile (either natural or synthetic) is used as the raw material. In a fluidized bed, the ore is chlorinated at ca. 900 °C under reducing conditions. The ore consists of 94-96% TiO2, the rest being oxides of other metals. The impurity chlorides from these other cations become separated in a series of condensers, and the pure TiCl4 is then evaporated, heated, and fed to a tubular or shaft burner together with a flow of superheated oxygen. As explained in Appendix D, a vigorous oxidation reaction occurs between the tetrachloride and oxygen, producing solid titania particles and gaseous chlorine. The temperature of the reaction mixture increases from ca. 900 °C to over 1400 °C at its maximum. In the hot reaction mixture, some gaseous intermediates might also be present. The activation energy of the reaction is quite high (ca. 88 kJ/mol at 1170 K), and thus, preheating of the reactant streams is essential. In an industrial burner, the ignition takes place at ca. 1200 K (mixed reactant temperature). Furthermore, the oxidation is an exothermic reaction with ∆H ) -177 kJ/mol (1273 K). The Gibbs energy change is -102 kJ/mol at 1273 K and -54 kJ/mol at 1773 K. Thus, an increase in the temperature of the reaction mixture might become a limiting factor for the overall conversion of titanium(IV) chloride to TiO2. To avoid excessive heating of the mixture, the reactor is provided with a cooling jacket with water circulation. Thermodynamic studies, e.g., that performed by Karlemo et al.,17 show that gaseous chlorine dissociates significantly at temperatures exceeding 1373 K. This dissociation is an endothermic process and, thus, is one additional factor affecting the temperature of the reaction mixture. The salient features of the burner chemistry in terms of the thermochemical model applied in the calculations are described in Appendix D. Once again, the model was constructed by combining the measured overall reaction kinetics with the thermodynamic multicomponent model, so as to produce a process model in which the conversion profiles of the reactor operating at a stationary production rate could be presented not only in terms of the temperature and concentrations but also in terms of the thermodynamic state functions for the calculated intermediate states. The conversion of the time scale of the reaction rates to a one-dimensional axial model in terms of the residence time was done using conventional steady-state assumptions for plug flow. The temperature of the reaction mixture and the temperature of the cooling water on each step are determined by applying a radial heat-transfer model in a temperature target calculation. The range of activation energies measured by Pratsinis and co-workers ranged from 73 to 101 kJ/mol18 when different O2/TiCl4 ratios were used in a laboratory experiment. The stoichiometric ratio O2/TiCl4 ) 1 gave a value of 73 kJ/mol, although a definite trend depending on the ratio of reactant inputs could not be established in their study. The kinetic data deduced from laboratory experiments must often be reevaluated against measured process data when a simulation model is being built. The earlier simulation data of Koukkari1 and Koukkari and Niemela¨19 were calculated using the average activation energy of 88 kJ/mol. In the recent modeling work performed by Koukkari, Keegel, and Penttila¨,20 different values of the activation energy Ea
2936
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002
maximum in the entropy production is visible, as the conversion from TiCl4 to TiO2 can reach 100% only at temperatures below 1573 K.
Figure 7. Entropy production and TiO2 conversion in the nonisothermal TiCl4 oxidation model.
were tried in the range of 71-102 kJ/mol. The best agreement with the measured heat-transfer and temperature data could actually be achieved with low Ea values, and thus, Ea ) 71 kJ/mol was adopted for the actual process model. The measured gas temperatures (published earlier19) were also in fair agreement with the modeled results. From the TiO2 yield curve (cf. appendix D), one can see that the reaction is slightly retarded as a result of the thermodynamic conversion limit at high temperatures. Additional heat transfer, which occurs in the cooling tube, completes the reaction. In practice, the location of the reaction exotherm can also be utilized for a secondary TiCl4 reactant injection. To evaluate the thermodynamic feasibility of the model, the entropy production was also calculated with different assumptions. The system is nonisothermal and contains no diffusive side streams (all of the material is conserved in the volume elements of the reactor tube, which are considered as closed systems during each step of simulation), and thus, one can see from eqs 8, 9, and 12 that diS can be obtained from dS by deducing the enthalpy-related heat-transfer contribution, dQ/T. Both the reaction rate model and the heat-transfer model have an effect on the entropy production curve. Figure 7 shows the thus-obtained entropy production for three values of the activation energy Ea. It can be seen that the rate of reaction combined with the heat-transfer model from each successive volume element must be such that positive entropy production is ensured for each step of the multistage model. After a slight decrease during the induction period of slow reaction rate, the entropy production increases until a maximum is reached, where the overall reaction rate also has a maximum and the slope of the conversion curve is steepest. In the cooling section, a secondary
Discussion The chemical change, as modeled with the thermochemical method in the examples, occurs via a series of successive “virtual” states. The system in each intermediate state is not in chemical equilibrium, as the chemical reactions leading toward its attainment have been suppressed by the calculation technique. However, the minimum of the Gibbs energy is calculated for each intermediate state, and for these intermediate states, one can use the phrase “frozen equilibria” as proposed by Guggenheim.21 The several chemical species present in the multicomponent system are virtually independent in the frozen states, and a chemical potential µ can be asssigned for each species. The thermochemical method provides a means to calculate the properties of such virtual states. The gradual attainment of the global equilibrium, represented by successive minima of the Gibbs energy, which are subdued to the measured reaction rate and other process conditions (heat and mass tranfer), also provides a simulated route for the occurring natural process. The combination of the reaction rates with the Gibbsian minimization is, to a degree, arbitrary, as exact information about the reaction mechanisms of the multicomponent reactive system is not often available. However, with the method presented above, one can easily follow the behavior of the thermodynamic functions during the chemical change and use the calculation as a guide to evaluate the physical consistency of the process model. In isothermal processes, this follows inherently, as the stepwise Gibbs energy minimization is directly deduced from the second law and thus leads to the monotonically descending Gibbs energy curves. These curves also define the entropy production of the isothermal process. For a nonisothermal system, the Gibbs energy change becomes defined by the applied reaction rate constraints and by the heat-transfer model, and the calculation provides a consistent entropy production profile throughout the reaction path. From the entropy production curve, the most active reaction zone with different model parameters and with different reactor structures can be deduced. For many cases, the thermodynamic state functions, which are directly obtained from the Gibbsian models, provide valuable information about the practical operation of the reactor. For example, in the recent work of Kjelstrup et al.,22 it was shown that the evaluation of the entropy production in chemical reactors can be a helpful tool in the optimization of a reactor process. The thermochemical method allows for the direct calculation of the entropy production in a multistage process model. Although the examples in this work were confined to the combination of an overall reaction rate process with the Gibbs energy method, more detailed methods to combine chain-reaction mechanisms with the multicomponent procedure can be applied. Similar thermodynamic conclusions should be valid for all such approaches. Appendix A. Lagrange Method for Gibbs Energy Minimization A number of methods can be used to minimize the Gibbs energy of multicomponent thermochemical systems. In the following, the method of Lagrange
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002 2937
undetermined multipliers is described for a simple system that consists of pure condensed phases and an ideal gas phase. The treatment is essentially that of Eriksson3 and has also been presented for systems with several mixture phases4 by the same author. The system consists of s condensed substances and of a gas phase with m constituents. The mole numbers and mole fractions of the constituents of the gas phase are denoted by ngk and xgk, respectively, and the constituents of the condensed phases by nck. The Gibbs energy of the system can be written as m
G)
∑
ngk
k)1
[( ) µ0
g
RT
k
+ ln p + ln
xgk
]
s
+
∑
k)1
nck
( ) µ0
c
RT
k
(A.1)
The number of elements in the system is l. The mass balance constraints for the elements of the system are m
The mass balance constraints then become defined by the single vector function ψ(n)
(∑
m+s
ψ(n) )
(j ) 1, ..., l)
n ) (n1, ..., nm+s) ) (ng1, ..., ngm, nc1, ..., ncs)
( ) ∂nk
-
T,p,ni
λjakj ) 0 ∑ j)1
( )
+ ln p + ln
( )
-
µ0
RT
k
∑
( ) µ0
RT
m
+ ln p
k
∑
m
nk +
k)1
∑
nk ln
k)1
( )
nk m
(A.3)
∑ nk
k)1
The mass balance constraints can, correspondingly, be presented as an (m + s) × l matrix A
ag11 . . . ag A ) m1 ac11 . . . g as1
. . . . . . . . . . . . . . . . . .
g a1l
agml g a1l
agsl
The total amount of each element is given in vector form as b ) (b1, ..., bl).
N
l
-
akjλj ) 0, ∑ j)1
k ) 1, ..., m
k
l
akjλj ) 0, ∑ j)1
k)1
nk
nk
(A.7)
µ0
RT
)
(A.6)
where λ1, ..., λl are the Lagrange undetermined multipliers. By performing the derivatives of the F function, the following set of equations is obtained
The Gibbs energy can then be written as a function of the molar amounts (n)
k)1
(k ) 1, 2, ..., m + s) (A.5)
ψ(n) ) 0
m+s
m+s
) ATn - b ) 0 (A.4)
l
∂F
0 0g 0c 0c µ0 ) (µ01, ..., µm+s ) ) (µ0g 1 , ..., µm , µ1 , ..., µs )
RT
l j)1
Minimization of the Gibbs energy is equivalent to minimization of the function F(n) with the limiting condition of ψ(n) ) 0. The respective Lagrange function is formed from eqs A3 and A4. The minimum of the Lagrange function is then obtained when all of its partial derivatives with respect to the molar amounts (nk) are zero. Concisely, the equations to be solved are given by
(A.2)
where aqkj is the number of atoms of the jth element in a molecule of the kth substance in the system (superscript q refers to the phase) and bj is the amount of the jth element in the system. The mass balances are the necessary constraints for the minimization problem, and eq A2 corresponds to the condition in eq 2 of the article text. The amounts of substances and the chemical potentials of the standard state can be expressed more concisely as
F(n) )
)
s
∑ agkj ngk + k)1 ∑ ackj nck ) bj k)1
G
akjnk - bj
k)1
k ) m + 1, ..., m + s
∑ akjnk - bj ) 0,
j ) 1, ..., l
(A.8)
(A.9)
The total number of moles in the gas phase is given m by N ) ∑k)1 nk. The number of equations is m + s + l, which include m + s + l variables, viz., nk with k ) 1, ..., m + s and λj with j ) 1, ..., l. Because of the thermochemical generality of the minimization problem, a wide variety of equilibrium systems can be treated with the said equations. However, even though eqs A.8 and A.9 are linear, there are no efficient algorithms that could be generally used for the nonlinear eqs A.7. The linearization of eqs A.7 can be done by applying Taylor series for the logarithmic terms, using the molar amounts nk ) nk(0) as the value points. These value points are used as initial guesses for the Taylor expression. With an iterative solution method, a series of calculations toward a final convergence point can then be performed. The calculations in this work were done with the ChemApp23 program. The matrix A defines the stoichiometric relations between the constituents of the system. This matrix can be extended to include constraints that are deduced from reaction kinetic restrictions.5,6 Such constraints might be due to the conservation of a metastable group, such as an aromatic ring in organic systems, but they might as well be due to the conservation of a constituent in the reactive processes. The modification of the matrix A brings about “kinetic control” of the number(s) of moles of the desired group(s) or substance(s) during the minimization procedure. The total mass of the multi-
2938
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002
component system is always conserved by applying the condition in eq 4 of the article text. Appendix B. Chemistry of the Multicomponent Redox System (See Table B1 for the stochiometry matrix of the Fe redox system.)
Overall reaction 2FeSO4 + 1/2O2 + H2SO4 T Fe2(SO4)3 + H2O (B.1) Redox reaction Fe2+ T Fe3+ + e-
(B.2)
Reaction rate (second-order, e.g., Iwai et al., 1982) 1 ) k′t + constant; k′ ) kpO2 mFe2+ Redox potential E ) E0 +
( )
aFe3+ RT ln zF aFe2+
(B.3)
(B.4)
aFez+ ) γFez+mFez+; γFez+ ) f(mk, T) The second-order reaction rate gives a straight line when 1/m(Fe2+) is plotted vs time. In this work, the concentrations were followed by taking samples, from which the Fe2+ concentrations were determined with the K2Cr2O7 method. The agreement of the measurements with the second-order trend were generally good, as indicated by Figure B1. The regression line log(k’) vs log(pO2) is also shown (Figure B2). When the pressure was varied between 5 and 40 bar (O2), the average value of the rate constant was obtained as k ) 19.5 × 10-7 mol-1 dm3 Pa-1 min-1 (353.15 K , 0.5 M H2SO4, 0.2 M FeSO4). The redox potential was followed with a Pt indicator coupled to a Ag/AgCl electrode through a KCl junction. The leading redox couple was assumed to be ferrous/ ferric ions. The thermodynamic multicomponent model was used to calculate the activities of the solute ions (ak) in terms of the extent of the redox reaction as given by the Nernst equation above (eq B.4) for Fe2+ and Fe3+ ions. The activity coefficients (γk) of the solute species (k) were calculated by using the Pitzer model for the specific ion-ion interaction parameters.14 The calcuTable B1. Stochiometry Matrix of the Fe Redox System N2(g) O2(g) H2(g) H2O(g) H2O H(+aq) OH(-aq) SO4(-2aq) HSO4(-aq) Fe(+2aq) Fe(+2aq)a Fe(+3aq) O2(aq) Fe(OH)2 Fe(OH)3 FeSO4 H2SO4
N
O
H
EA
S
Fe
Fea
2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 2 0 1 1 0 1 4 4 0 0 0 2 2 3 4 4
0 0 2 2 2 1 1 0 1 0 0 0 0 2 3 0 2
0 0 0 0 0 -1 1 2 1 -2 -2 -3 0 0 0 0 0
0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 1 1
0 0 0 0 0 0 0 0 0 1 0 1 0 1 1 1 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0
a Reaction rate constraint applied for Fe2+ oxidation in the Gibbs energy model.
Figure B1. Change of Fe(II) concentration with oxygen pressure. Plot of 1/m(Fe2+) vs time for the reaction 2FeSO4 + 1/2O2 + H2SO4 T Fe2(SO4)3 + H2O.
Figure B2. Plot of log(k’) vs log(pO2) for the reaction 2FeSO4 + 1/ O + H SO T Fe (SO ) + H O. 2 2 2 4 2 4 3 2
lated potential was converted to the SHE (standard hydrogen) scale by using the empirical correction for the KCl junction given by MacDonald et al.12 At 353.15 K (80 ˚C), ∆ECORR ≈ 0.196 V when the KCl concentration is 0.5 M. Appendix C. The Dissolution of CaCO3 by Sulfuric Acid When sulfuric acid is added to an aqueous solution containing undissolved calcium carbonate (Figure C1), the overall reaction produces soluble calcium sulfate and carbon dioxide through the following sequence
H2SO4 f H+ + HSO4-
(C.1)
CaCO3 + H+ + HSO4- f Ca2+ + HCO3- + HSO4(C.2) HCO3- + HSO4- f H2O + SO42- + CO2v
(C.3)
CaCO3 + H2SO4 f Ca2+ + SO42- + H2O + CO2v Subsequent acidifications in the solution lead to successive pH drops, which are then offset by the presence
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002 2939
and the cooling water are combined into a steady-state simulation in terms of the reactor axis (x) (see Figure D1). The temperature of the reaction mixture (T) and the heat recovered in different axial sections of the reactor tube can be measured.19,20 The key reactions and their data (at 1273 K) are as follows:
Figure C1. Schematic of an aqueous solution containing undissolved calcium carbonate.
TiCl4(g) + O2(g) T TiO2(s) + 2Cl2(g)
(D.1)
Cl2(g) T 2Cl(g)
(D.2)
Reaction
∆H/kJ mol-1
∆G/kJ mol-1
Ea/kJ mol-1
D.1 D.2
-177 250
-102 98
71-102 -
The reaction rate parameters for the consumption of TiCl4 in the oxidation reaction were determined by Pratsinis et al. in 1989.18 The overall rate equation was given as
d[TiCl4] ) -(k + k′x[O2])[TiCl4] dt
(D.3)
with
k ) A e-Ea/RT, k′ ) A′ e-Ea/RT
Figure C2. Effect on pH of the addition of 1 M sulfuric acid to a stirred solution of tap water and 2.49 g of CaCO3.
of dissolving CaCO3. Figure C2 shows the effect on pH of three additions (8 mL each) of 1 M sulfuric acid to a stirred solution of tap water (5000 mL) to which 2.49 g of CaCO3 was added in the beginning of the experiment. The Ca2+ content of the solution was simultaneously followed with a Ca2+ ion-selective electrode15 (Orion M 97/20). Appendix D. One-Dimensional Model of the TiCl4 Burner The extent of reaction (ξr), the enthalpies of reaction, and the heat exchange between the reaction mixture
Figure D1. One-dimensional model of the TiCl4 burner.
where k and k′ are given in the usual Arrhenius form with the numerical values for the frequency factors A and A′ and the activation energy Ea as deduced by Pratsinis et al. The values they determined are A ) 8.26 × 104 s-1, A′ ) 1.4 × 105 (dm3 mol-1)1/2 s-1, and Ea ) 73 kJ mol-1. One can conclude that, at close-to-stoichiometric compositions, k′[O2]1/2 , k, and the reaction rate is practically first-order in terms of the TiCl4 concentration. In addition to chlorine association-dissociation reactions, other side reactions such as oxychloride formation take place in the hot regions of the burner. Although these reaction processes have a significant effect on the total enthalpy change, their reaction rates are not known. Thus, Cl2 dissociation and other side reactions were assumed to be fast (in equilibrium), and the reactor model was constructed by combining the global reaction rate with the multicomponent Gibbsian calculation of a TiCl4-O2 system.1,5 The temperature of each successive reactor volume element was then obtained with an iterative target calculation from estimated heat-transfer
2940
Ind. Eng. Chem. Res., Vol. 41, No. 12, 2002
Figure D2. Sectional heat losses for predefined segments along the reactor’s longitudinal axis.
Figure D3. Temperature profiles of various components of the reactor.
data. The heat transfer from the hot reaction mixture to the surrounding cooling water can be described as a combination of radiation and convection in the radial direction.5,20 The typical profiles calculated from the reactor model include the sectional heat losses for predefined segments along the reactor’s longitudinal axis (Figure D2) and the temperature profiles of various components of the reactor (Figure D3). Literature Cited (1) Koukkari, P. A physicochemical method to calculate timedependent reaction mixtures. Comput. Chem. Eng. 1993, 17 (12), 1157-1165. (2) Eriksson, G. Thermodynamic studies of high-temperature equilibria III. Acta Chem. Scand. 1971, 25 (7), 2651-2656. (3) Eriksson, G. Thermodynamic studies of high-temperature equilibria. XII. SOLGASMIX, a computer program for calculating equilibrium compositions in multiphase systems. Chem. Scr. 1975, 8, 100-103. (4) Smith, W. R.; Missen, R. W. Chemical Reaction Equilibrium Analysis: Theory and Algorithms; Krieger Publishing Company: Malabar, FL, 1991. (5) Koukkari, P. A physicochemical reactor calculation by successive stationary states. Acta Polytech. Scand., Chem. Technol. Ser. 1995, 224. (6) Koukkari, P.; Pajarre, R.; Hack, K. Setting kinetic controls for complex equilibrium calculations. Z. Metallkd. 2001, 92 (10), 1151-1157.
(7) Koukkari, P.; Laukkanen, I.; Liukkonen, S. Combination of overall reaction rate with Gibbs energy minimization. Fluid Phase Equilib. 1995, 136, 345-362. (8) Krupp, U.; Christ, H.-J. Internal nitridation of nickel-based alloys. Part II. Oxid. Met. 1999, 52 (3/4), 299-320. (9) Robertson, D. G. C. The computation of the kinetics of the reactions between multiple phases. In TMS EPD Congress; Warren, G. W., Ed.; TMS: Warrendale, PA, 1995; pp 347-361. (10) Kondepudi, D.; Prigogine I. Modern Thermodynamics: From Heat Engines to Dissipative Structures; John Wiley & Sons: New York, 1998. (11) Pohjola, K. The oxidation-reduction potential of Fe in pressurized conditions. M.Sc Thesis, Department of Materials Technology and Mining Technology, Helsinki University of Technology, Espoo, Finland, 1997 (in Finnish). (12) McDonald, N.; Scott, A.; Wentreck, P. External reference electrodes for use in high-temperature aqueous systems. J. Electrochem. Soc. 1979, 126 (6), 908-911. (13) Iwai, M.; Majima, H.; Awakura, Y. Oxidation of Fe(II) in sulfuric acid solutions with dissolved molecular oxygen. Met. Trans. 1982, 13B, 311-318. (14) Pitzer, K. S. Thermodynamics, 3rd ed., McGraw-Hill; New York, 1995. (15) Koukkari, P.; Pajarre, R.; Pakarinen, H.; Salminen, J. Practical multiphase models for aqueous process solutions. Ind. Eng. Chem. Res. 2001, 40, 5014-5020. (16) Compton, R.; Daly, P. The dissolution/precipitation kinetics of calcium carbonate: An assessment of various kinetic equations using a rotating disk method. J. Colloid Interface Sci. 1987, 115 (2), 493-498. (17) Karlemo, B.; Koukkari, P.; Paloniemi, J. A Kinetic and Thermodynamic Study of the Oxidation of Titanium(IV) Chloride at High Temperatures; Report TKK-V-B70; Helsinki University of Technology: Espoo, Finland, 1992. (18) Pratsinis, S. E.; Bai, H.; Frenkloch, M.; Mastrangelo, S. V. R. Kinetics of titanium(IV) chloride oxidation. J. Am. Ceram. Soc. 1990, 73 (7), 2158-62. (19) Koukkari, P.; Niemela¨, J. Time-dependent reactor simulation by stationary state thermochemistry. Comput. Chem. Eng. 1997, 21 (3), 245-253. (20) Koukkari, P.; Keegel, M.; Penttila¨, K. Coupled thermodynamic and kinetic models for high-temperature processes. In HighTemperature Materials Chemistry; Hilpert, K., Froben F. W., Singheiser L., Eds; Schriften des Forschungszentrum Julich; IUPAC: Research Triangle Park, NC, 2000; Vol. 15, Part I, pp 253-246. (21) Guggenheim, E. A. Thermodynamics, 6th ed.; NorthHolland Publishing Company: Amsterdam, 1975. (22) Kjelstrup, S.; Johannessen, E.; Ro¨sjorde, A.; Nummedal, L.; Bedeaux, D. Minimizing the entropy production for the methanol producing reaction in a methanol reactor. Int. J. Appl. Thermodyn. 2000, 3 (4), 147-153. (23) Eriksson, G.; Hack, K.; Petersen, S. ChemAppsA programmable thermodynamic calculation interface. In Werkstoffswoche ‘96, Symposium 8: Simulation Modellierung, Informationsysteme; DGM Informationsgesellschaft mbH: Frankfurt, Germany, 1997; p 47.
Received for review June 6, 2001 Revised manuscript received February 1, 2002 Accepted February 6, 2002 IE010498T
