A reliable procedure to predict salt precipitation in pure phases

THERMODYNAMICS AND SEPARATION PROCESSES

A reliable procedure to predict salt precipitation in pure phases

P. S. O. Beltrão^I; M. J. E. M. Cardoso^II; M. Castier^III,*

^IBuckman Laboratories, Campinas - SP, Brazil

^IIInstituto de Química, Universidade Federal do Rio de Janeiro, Rio de Janeiro - RJ, Brazil

^IIIDepartment of Chemical and Petroleum Engineering, U.A.E. University, Al-Ain, P.O. Box: 17555, United Arab Emirates. E-mail: [e-mail]

ABSTRACT

This article proposes a new procedure to compute solid-liquid equilibrium in electrolyte systems that may form pure solid phases at a given temperature, pressure, and global composition. The procedure combines three sub-procedures: phase stability test, minimization of the Gibbs free energy with a stoichiometric formulation of the salt-forming reactions to compute phase splitting, and a phase elimination test. After the phase splitting calculation for a system configuration that has a certain number of phases, the phase stability test establishes whether including an additional phase will reduce the Gibbs free energy further. The criteria used for phase stability may lead, in some cases, to the premature inclusion of phases that should be absent from the final solution but, if this happens, the phase elimination sub-procedure removes them. It is possible to use the procedure with several excess Gibbs free energy models for liquid phase behavior. The procedure has proven to be reliable and fast and the results are in good agreement with literature data.

Keywords: Thermodynamics; Equilibrium; Aqueous solutions; Solubility; Precipitation.

INTRODUCTION

Knowledge of the thermodynamic properties of electrolyte solutions is relevant for the design of many classes of chemical processes (Chen, 2006), some of which of very large scale. Freshwater, for example, comes mostly from desalination plants in some countries because of its scarcity from natural sources. Processes in the oil and gas industry are also generally of large scale and the properties of electrolyte systems are relevant in certain situations as in sour gas sequestration problems (Ferrando et al., 2006) and oil-brine separations. The food industry is another example of a sector that needs to design and operate processes in which electrolyte solutions are present (Bruin, 1999).

This article focuses on salt solubilities. Solid formation is sometimes desired, as in crystallization processes, but is a nuisance in other situations, as when mineral scales form inside chemical processing equipment. Phase diagrams that present solubility data are available in the literature for many systems. Direct use of experimental data or model predictions (Chen, 2006; Chiavone and Rasmussen, 2000; Lin et al., 2007; Messnaouia et al., 2008), as in the IVC-ELEC program for phase diagrams of electrolyte systems (Thomsen, 1997), allow their preparation. These diagrams provide insight about the complex phase behavior some systems exhibit and are useful to set proper operation ranges for process variables. For process design, it is also of interest to find the equilibrium conditions of an electrolyte system of known global composition, temperature, and pressure. Recent developments in this subject include the work of Montastruc et al. (2004) who studied the precipitation of calcium phosphate at different pH and formulated the problem as the minimization of the Gibbs free energy of the system, which is carried out by a genetic algorithm in the initial iterations, followed by a successive quadratic programming method in the final ones. Cezac (2005) formulated and solved the salt precipitation problem using a set of non-linear equations, which explicitly included the elemental balance and electroneutrality equations, i.e., adopting a non-stoichiometric approach.

Problems with specifications of temperature, pressure, and global mole fractions, when referring to fluid phase separations, are known as isothermal flash problems. The equilibrium conditions in an isothermal flash problem minimize the Gibbs free energy of the system. In some specialized situations the number of phases is known beforehand, but in general situations, an isothermal flash procedure has to find the number of phases that coexist and the number of moles of each component in each phase. An approach is to overestimate the number of phases present and let the component mole numbers in absent phases approach zero as the numerical calculations proceeed. The disadvantage is that the phase equilibrium problem is set up with a number of unknowns that is, in general, larger than it should be. Another approach is to start computations with a single phase and add phases, one by one, as dictated by the outcome of the global phase stability test (Michelsen, 1982a). Experience with this second approach (Michelsen 1982a,b; Castier et al., 1989) indicates that the global phase stability test may include phases prematurely, which are absent from the equilibrium state. Therefore, a general flash procedure that adds phases based on the global phase stability test needs to have a procedure to remove phases included prematurely (Michelsen 1982a,b; Castier et al., 1989). This is necessary not only in isothermal flashes but also in other types of calculations, such as isochoric-isoenergetic flash computations (Castier, 2009).

A numerical technique to solve the isothermal flash problem is to use a minimization algorithm. Another is to solve the set of non-linear equations that results from one of the many possible ways of rewriting the equilibrium conditions. It is possible to extend any of these formulations to account for chemical reactions by using either the non-stoichiometric or the stoichiometric approach. The former approach imposes mass conservation via elemental balances (and electroneutrality equations in the case of electrolyte systems) explicitly while the latter uses component balances via stoichiometric coefficients and extents of reaction instead. Different flash algorithms result from the choice of numerical technique and approach to incorporating chemical reactions, among other details.

It is, in principle, possible to adapt these procedures to solve isothermal flash-like problems involving electrolyte systems and there are several examples of vapor-liquid equilibrium computations for electrolyte systems that use Gibbs free energy minimization (Cherif et al, 2000, 2002). An example with solid phases is a Gibbs free energy minimization method that uses simple thermodynamic phase models (ideal K-values) to predict gas hydrate precipitation (Ballard and Sloan, 2004) in systems that contain salts. Other examples include the SOLGASMIX program, which minimizes the Gibbs free energy and uses the Pitzer model for electrolyte solutions, and the ESP program, which solves the equilibrium algebraic equations and uses the Bromley model for electrolyte solutions. The Oak Ridge National Laboratory compared these two programs (Beahm et al., 2000; Toghiani et al., 2000) focusing on the effect of formation properties and thermodynamic models on results rather than on numerical performance. SOLGASMIX later became part of the FactSage package (Bale et al., 2002) of thermochemical software and databases. Educational software with graphical interface for simple salt solubility problems is also available (University of Colorado, 2009; Gil and Paiva, 2006). A formulation based on solving the algebraic equations that relate the solubility products to the activities has been recently used to predict the formation of sulfate and carbonate scaling minerals (Villafáfila García, 2005, 2006), but details about the numerical implementation are scarce in these references.

Geochemists often perform chemical and phase equilibrium calculations in dilute aqueous systems and refer to these computations as speciation in spite of IUPAC's criticism (Templeton et al., 2000) of the term. Lothenbach et al.(2007) mentions the following geochemical codes for equilibrium computations: MINEQL++ (2009) and CHESS (van der Lee, 1998; van der Lee and DeWind, 2002), which are commercial software, and GEMS (GEM-Selektor) and MINTEQA2, which are freely available, at least for academic use. Kulik (2007) describes the GEM-Selektor package, which is a large set of tools for equilibrium calculations for geochemical applications that computes activity coefficients by using models such as those of Debye-Hückel, Pitzer, and Brønsted-Guggenheim. GEM finds the equilibrium conditions via minimization of thermodynamic potentials under different specifications using a constrained optimization method. MINTEQA2 (Allison et al., 1991) is a geochemical software developed by the United States of America Environmental Protection Agency (EPA). MINTEQA2 uses the Debye-Hückel or the Davies formulas for activity coefficients and EPA does not recommend using the program at ionic strengths higher than 0.5. MINTEQA2 finds the equilibrium conditions by solving the set of equilibrium equations. CHESS (van der Lee, 1998; van der Lee and DeWind, 2002) uses the Debye-Hückel or different version of the Davies formulas for activity coefficients and solves the set of equations that represent the equilibrium conditions. It uses a modified form of the Newton-Raphson method intented to increase the probability of achieving numerical convergence (van der Lee and DeWind, 2002). CHESS adds one solid phase at a time, based on a saturation index, and removes solid phases when they become undersaturated. MODELm (Huber et al., 2002) is another example of code to compute the equilibrium condition of aqueous systems with trace amounts of ions and organic compounds. The numerical technique is to solve the set of equations that represents the equilibrium conditions.

This paper presents a new procedure to compute SLE in aqueous electrolyte systems of specified temperature, pressure, and initial amounts of water, cations, and anions. Equivalent to specifying ion amounts, it is possible to specify the amounts of salts added, as these quantities can be easily converted to equivalent ion-based specifications. For example, it is equivalent to specify an initial amount of 1 mol of Na₂SO₄ or 2 moles of Na⁺ and 1 mol of . Given these specifications, the procedure finds which species will actually form, the number of phases present, and how the formed species distribute among these phases (i.e., performs a speciation calculation).

The equilibrium problem is formulated and solved as the minimization of the Gibbs free energy of the system. The solution procedure takes advantage of developments in the area of flash calculations for non-elecrolyte systems. Unlike most of the geochemical codes reviewed here, calculations are not limited to small values of ionic strength because of the use of the extended UNIQUAC model. However, the procedure is general and can use several activity coefficient models, including those used in geochemical software. Unlike the work of Montastruc et al. (2004), we do not use an empirical optimization method, such as the genetic algorithm. The type of problem we solve, combined with the structure of our algorithm, guarantees that our procedure finds the global minimum of the Gibbs free energy by only using a local optimization algorithm. Differently from the work of Cezac (2005), we use the stoichiometric approach because it satisfies the elemental balances and electroneutrality equation automatically with a smaller number of unknowns to be solved by the numerical method. Furthermore, our formulation allows the use of an unconstrained optimization to minimize the Gibbs free energy, in contrast to the use of constrained optimization methods, generally more complicated, as done in some of the references reviewed here. In addition, we use the global phase stability test to decide about the addition of solid phases, one at a time, thereby finding the correct number of phases in the system while avoiding the need to start calculations with a large number of them. A potential risk of this approach, especially close to phase boundaries, is that the stability test may include a solid phase that should be absent from the problem solution. To avoid this, during the phase splitting part of the computation, we try to merge the solid phase present with the smallest number of moles with the liquid phase. We accept this action if it reduces the Gibbs free energy of the system and then decrease the number of phases by one. The combination of steps for automatic phase addition and removal and of a fast numerical method to minimize the Gibbs free energy in the phase splitting part of the algorithm gives rise to a new and robust procedure for the solution of SLE problems whose details and applications are discussed in the next sections.

PHASE EQUILIBRIUM CALCULATIONS

The major assumptions are:

1. the liquid phase contains solvent(s) and ions but does not contain undissociated salts;

2. each solid phase contains a single substance, which is an undissociated salt, hydrated or not.

The overall computation procedure is as follows. We initially assume that a single liquid phase contains all the solvent(s) and ions. We then test the stability of this liquid phase with respect to the possible formation of solid phases. If it is unstable, we add one new pure solid phase and a new chemical reaction to obtain the molecular form of the precipitating salt from its ions. The next step is to minimize the Gibbs free energy of the system to find the amount of each species in each phase (phase splitting). We repeat these cycles of phase stability and phase splitting calculations until the system is found to be stable. We also test for the possible elimination of solid phases during the phase splitting step.

The next subsections discuss details of the thermodynamic description used for each phase and of our phase equilibrium calculations, including the problem formulation and numerical procedure. This section finishes with a subsection that discusses the ability of the algorithm to find the global minimum of the Gibbs free energy.

Thermodynamic Description of the Phases

We assume that precipitated salts form pure phases. The procedure is, however, general with respect to the liquid phase and can use activity coefficients derived from many excess Gibbs free energy expressions. In this work, we use the extended UNIQUAC model for electrolyte solutions. There are slightly different versions and parametrizations of this model in the literature (Nicolaisen et al., 1993; Pinho and Macedo, 1996; Thomsen et al., 1996). We used the model version and the parameters reported in the work of Nicolaisen et al. (1993), along with the necessary Gibbs free energies and enthalpies of formation and heat capacities at constant pressure for the ions, salt, and water from the same source (Nicolaisen et al., 1993).

Chemical potentials are the most important properties for the computations discussed in this article. For a salt S in its pure solid phase, the chemical potential (µ_S) is:

where R and T are the universal gas constant and the absolute temperature, respectively, is the molar Gibbs free energy of formation of the pure solid salt at the system temperature, and is the salt activity in the pure solid phase, which is equal to 1. We therefore have that µ_S= .

For water in the liquid phase, the chemical potential () is:

where is the molar Gibbs free energy of formation of water as a pure liquid at the system temperature and , , and are the activity, mole fraction, and activity coefficient of water in the liquid phase, respectively. Extended UNIQUAC gives activity coefficients on a mole fraction scale and uses the symmetric convention for solvents. Therefore, values of calculated by extended UNIQUAC can be used directly in Eq. 2.

For an ion I in the liquid phase, the chemical potential is:

where is the molar Gibbs free energy of formation of ion I compatible with the use of unsymmetrically normalized activity coefficients (denoted as ) based on the mole fraction scale. This is the type of activity coefficient extended UNIQUAC computes for ions, which can therefore be directly used in Eq. 3.

The Gibbs free energies and enthalpies of formation and heat capacities at constant pressure reported by Nicolaisen et al. (1993) allow direct computations of and , as will be detailed later. For ions, however, they used:

where m_I is the molality of ion I and is the molar Gibbs free energy of formation of ion I compatible with the use of unsymmetrically normalized activity coefficients () based on the molal scale.

Comparing Eqs. 3 and 4, it results that:

The left-hand side of Eq. 5 does not depend on the mole fraction or molality of ion I in the system because is a formation property whose value depends on the formation state and not on the state of the system. Therefore, the right hand side (RHS) must depend on the formation state only. A convenient way to find the value of the RHS is to consider an aqueous solution whose amounts of solute tend to zero. In the mathematical formulas that follow, this will be denoted by solution → water. At this condition, the ratio m_I/x_I becomes:

where the summation index j runs over all components, ionic or not, except water, and is the molar mass of water, in . Furthermore, as the system becomes dilute, m_I→ 0 and x_I→ 0, such that → 1 and → 1 . As both values tend to 1, their ratio, /, also tends to 1.

From these considerations, Eq. 5 becomes:

The final issue related to the thermodynamic description of the phases is how to compute the values of , , and that appear in Eqs. 1, 2, and 7, respectively, using the data provided by Nicolaisen et al. (1993). Adopting as a general symbol to represent , , , or and under the assumption that the heat capacity at constant pressure does not depend on temperature, then:

where , , and are the values tabulated by Nicolaisen et al. (1993), and is the reference temperature, equal to 298.15 K.

Phase Stability

The tangent plane criterion establishes whether a phase of given composition, temperature (T) and pressure (P) (in which the chemical potential of each species i is µ_i) is stable or not. Phase stability is ensured when:

for every set of mole fractions in the trial phase, , in which the chemical potential of each species i is . Practical implementation of the test in general situations requires the location of the global minimum of the function Ω_G. If the global minimum value of Ω_G is not positive, the original phase is unstable, and the situation in which it is equal to zero corresponds to the stability limit.

In this work, we test the stability of the liquid phase with respect to the possible formation of solid phases that only contain pure undissociated salts. In this case, Inequality 9 simplifies to:

where subscript S refers to an undissociated salt. The chemical potential of the salt in the trial solid phase, , is:

Salt activity in the pure solid phase () is equal to 1 and we therefore have that:

where is the molar Gibbs free energy of formation of the salt.

The chemical potential of the salt in solution, µ_S, can be determined from the following general reaction:

In this equation, is the number of ionic species, (ION)_i is the chemical symbol of an ion of type , represents the juxtaposition of chemical element symbols in the salt formula. In a hydrated salt, ≠ 0, otherwise = 0. The λ coefficients of and of each ion i relate to the respective stoichiometric coefficients with sign convention (negative for reactants and positive for products), ν, as follows:

The chemical equilibrium condition is given by:

or equivalently by:

Substituting Eqs. 12 and 17 into Inequality 10 and writing the chemical potentials in terms of activities, we obtain the following form of the phase stability criterion:

If many types of cations and anions are present, they may form several possible salts, including different hydrated forms. In our implementation, it is necessary to provide input data for all salts that may precipitate in the system. Computations are sequential in the sense that, at each step, we test phase stability with respect to the possible formation of one additional solid phase. Instability of the original system triggers the phase splitting calculations. This alternating sequence of phase stability and phase splitting calculations continues until the system is found to be stable.

At each step, the value of the function Ω_G for each salt, calculated with the current values of the system properties, defines the order of attempted phase additions. If the system is unstable, the new phase added to the system is that with the smallest value of Ω_G (provided that Ω_G < 0), i.e., the new phase with the largest violation of the phase stability condition according to the Ω_G function. Thomsen (1997) used a solubility index, defined as the activity product of a given salt divided by its solubility product, as the criterion to rank candidate salts for precipitation, but Inequality 18 has a more direct connection with the primary phase equilibrium criterion, which is the minimization of the Gibbs free energy.

Phase Splitting

We assume that salts in molecular form are absent from the liquid phase and that, if formed, precipitate as pure solid phases. The SLE of a saturated aqueous solution with precipitated salts in pure phases can be formulated as a chemical equilibrium problem in which the type of reaction shown in Eq. 13 takes place.

Phase splitting is determined by minimization of the Gibbs free energy (G) of the system which, for convenience, we divide by the product RT. In a multiphase system with n_c species and n_p phases, we have that:

In the equations that follow, we use a general notation but simplications apply, as described here, to a system that has a single liquid phase and multiple solid phases, each of them containing a single precipitated salt. The number of moles of species i in phase j is:

where θ_ij represents the fraction of species i present in phase j. In these calculations, only the extents of the reactions that form the precipitating salts are iterated variables because the distribution of species among the phases is known a priori: each salt formed in the system will be present only in its corresponding pure solid phase; ions and free water (not present in hydrated salts) will be present only in the liquid phase. Therefore, for the molecular form of a salt, θ is equal to one in its pure solid phase and zero elsewhere in the system. For the solvent(s) and ions, θ is equal to one in the liquid phase and zero elsewhere in the system. The symbol n_i represents the number of moles of species i in the system, given by:

where n_i0 is the initial number of moles of species i, n_r is the number of chemical reactions, v_ik is the stoichiometric coefficient of species i in reaction k, and ξ_k is the extent of reaction k. These extents of reaction are the independent variables in the minimization of the Gibbs free energy. Expressions for the first and second derivatives of the Gibbs free energy with respect to extents of reactions are available in the work of Castier et al. (1989). The terms of the gradient are:

The elements of the Hessian matrix are:

where a_ij is the activity of species i in phase j. Expressions for the derivatives of the activities with respect to number of moles are available in the work of Cardoso (1989).

We use the Murray local optimization method, described in a book (Murray, 1972) by the same author, to minimize the Gibbs free energy of the system with initial estimates for the extents of reaction generated as follows: in a partition that has (n_p + 1) phases, we use the converged values of the extents of reactions from the phase splitting computations with n_p phases as initial estimates for the pre-existing salts; for the new precipitating salt added by the phase stability test, we use an initial estimate equal to 10^-6 moles.

Phase Elimination

The sequential application of the stability test and phase splitting procedures may introduce a phase in an early stage of the calculation that needs to be removed later. A typical situation is the early inclusion of a phase whose amount gradually decreases as computations proceed and more solid phases are added. Phase elimination is not attempted for system partitions that have two phases, as this would take the system back to the single phase configuration already detected as unstable. Phase elimination is, therefore, attempted only for system partitions that have three or more phases.

Phase elimination proceeds as follows: at each iteration, we incorporate the solid salt present with the smallest number of moles (in the form of its originating ions) into the liquid phase and test whether this action decreases the total Gibbs free energy of the system. If it does, the corresponding solid phase is eliminated, otherwise the solid phase is kept. Furthermore, to reduce the risk of an undue phase elimination, we only perform phase elimination attempts after two iterations of the phase splitting procedure for a given number of phases (when n_p> 3).

Global Minimum of the Gibbs Free Energy

The use of globally convergent methods, either stochastic or deterministic, to solve phase and/or chemical equilibrium problems and related parameter fitting problems has increased in recent years. Commonly used stochastic methods include simulated annealing (e.g., Bonilla-Petriciolet et al., 2008), genetic algorithms (e.g., Alvarez et al., 2008), and particle swarm (e.g., Rahman et al., 2009). Deterministic methods generally include a procedure to subdivide the root search area and systematically discard regions that do not contain a root (e.g., Corazza et al., 2007, Simoni et al., 2007).

Many types of phase equilibrium algorithms exist, among which those that alternate phase stability and phase splitting calculations, as done here. If an unstable configuration is detected in the stability test, an additional phase is included in the system. With the updated number of phases, minimization of the Gibbs free energy provides results for phase splitting, i.e., phase amounts and compositions.

Detection of phase instability means that it is possible to reduce the Gibbs free energy by inclusion of an additional phase. Using adequate initial estimates and a local minimization algorithm that never allows increases of the Gibbs free energy from one iteration to the next, the minimization procedure lowers the Gibbs free energy in all iterations and stops at a local minimum.

The key point is to decide whether the local minimum of the Gibbs free energy is also its global minimum. The phase stability test should be reliably implemented, in such a way that, if the condition is not that of the global minimum, it will propose the addition of a new phase that will contribute to lower the Gibbs free energy further. Fully reliable implementations of the phase stability test are based on global optimization algorithms (Floudas and Gounaris, 2009), such as interval analysis, to minimize the function Ω_G (Eq. 9). Stochastic optimization methods can also be used (e.g., Nagatani et al., 2008), but they do not provide a mathematical guarantee of locating the global minimum of the phase stability function. In fact, if the phase stability test is carried out properly and always finds the global minimum of the function Ω_G, a local optimization algorithm can be used to minimize the Gibbs free energy during the phase split calculation. The combination of these two procedures guarantees that the global minimum of the Gibbs free energy is found because, if the local minimization method used for phase splitting does not locate the global minimum of the Gibbs free energy, the global optimization method used for the phase stability test detects this fact and includes a new phase.

The assumption of pure solid phases introduces considerable simplification because phase stability or instability can be reliably determined by just examining the signal of Inequality 10, without the need for any numerical method. It is only necessary to test for the possible formation of all salts that may form from the ions present in the system, in hydrated form or not.

In summary, even though the method proposed here only uses a local optimization algorithm, it finds the global minimum of the Gibbs free energy and, therefore, provides correct solutions to SLE problems containing pure solid phases.

RESULTS

We discuss examples with a single salt in water and with multiple cations and anions that cause the precipitation of several pure phases. In all examples, initial values for the unknowns were obtained according to the procedure described in the paper and no ad-hoc value was used. The convergence criterion is an Euclidian norm of the objective function gradient smaller than 10^-8. The primary result of our calculations is the number of moles of each species, solute or solvent, in each phase. For easier comparison with literature results, some of these primary results were converted to concentration scales upon completion of the equilibrium calculations. In particular, because the same excess Gibbs free energy model and parameters are used, it is possible to compare our results directly with those of Nicolaisen et al. (1993).

Precipitation of Pure Salts

We computed the effect of temperature on the solubility of different pure salts in water. Because the procedure does not aim at calculating saturation points directly but is a flash-like procedure, it needs the specification of the global amounts of water and ions present. In these calculations, we specified an excess of ions compared to the amount expected to be soluble. As discussed in the previous section, the algorithm starts by assuming that water and the ions are all in the liquid phase and tests if this phase is stable. As these configurations were detected as unstable, the procedure computes the amount of salt precipitated and, therefore, the amount of ions left in the liquid phase corresponds to the salt solubility.

Tables 1 and 2 show solubilities of NaCl and K₂SO₄, respectively, in water as a function of temperature. The average absolute deviations between calculated and experimental values are equal to 1.6% for NaCl and 2.7% for K₂SO₄. We made analogous calculations for NaCl and K₂SO₄, not shown here. For the four salts, the results are in excellent agreement with those of Nicolaisen et al. (1993). These tables also include the activity coefficients based on the mole fraction scale for water (symmetric convention) and for the ions (unsymmetric convention) at each equilibrium condition. Computations of the SL partitioning typically converge in three iterations of the Murray method.

System Na₂SO₄ - K₂SO₄ - NaK₃(SO₄)₂ - Na₂SO₄.10H₂O - H₂O

The salts that precipitate in systems that have the Na⁺, K⁺, and ions in water depend on temperature. The experimental evidence is that systems at 323.15 K may form not only the simple salts ( Na₂SO₄ and K₂SO₄) but also the complex double salt NaK₃(SO₄)₂. At 298.15 K the salts that can precipitate are Na₂SO₄.10H₂O, K₂SO₄, and NaK₃(SO₄)₂. In the formulation of all cases reported here, we assume the possibility of forming the following four salts: Na₂SO₄, K₂SO₄, NaK₃(SO₄)₂, and Na₂SO₄.10H₂O. The procedure determines which of them, if any, precipitate(s) for a given set of case specifications.

Temperature of 323.15 K

Table 3 summarizes the specifications and results of several test cases at the temperature of 323.15 K. In all of them, the specified amount of water is 55.50844 moles (1 kg). Cases 1-3 illustrate situations in which a single salt precipitates: K₂SO₄, NaK₃(SO₄)₂, and Na₂SO₄, respectively.

Case 4 illustrates a more complex situation in which two salts precipitate: Na₂SO₄ and NaK₃(SO₄)₂. During the course of calculations, pure NaK₃(SO₄)₂ is the first solid phase added to the system, even though its number of moles is smaller than that of the pure solid phase of Na₂SO₄ in the final case solution, as shown in Table 3. In Case 5, two other salts precipitate: K₂SO₄ and NaK₃(SO₄)₂. As in Case 4, NaK₃(SO₄)₂ is the first pure solid phase added to the system during the calculations.

The phase diagrams of Nicolaisen et al. (1993) for this system report solubilities expressed as molalities of Na₂SO₄ and K₂SO₄. Once converted to this concentration scale, all the results of Table 3 are in excellent qualitative and quantitative agreement with those of Nicolaisen et al. (1993). Furthermore, the precipitation of Na₂SO₄.10H₂O is not predicted in any of these cases, in agreement with experimental evidence.

Phase splitting calculations typically required three or four iterations of the Murray method.

Temperature of 298.15 K

Table 4 presents specifications and results at 298.15 K. In all of them, the specified initial amount of water is 55.50844 moles (1 kg), but it should be observed that Na₂SO₄.10H₂O precipitates in Cases 7 and 9. When it does, the amount of water in the liquid phase decreases because part of the water goes to the solid phase. In Case 6, only K₂SO₄ precipitates. In Case 7, we swap the amounts of Na⁺ and K⁺ with respect to the values in Case 6, and the procedure predicts the precipitation of Na₂SO₄.10H₂O. With the specifications of Case 8, only NaK₃(SO₄)₂ precipitates. Cases 9 and 10 show the simultaneous precipitation of two salts. Once our results are converted to the molality scale used by Nicolaisen et al. (1993), agreement with their calculations is excellent. In agreement with experimental evidence, the procedure does not predict the precipitation of Na₂SO₄ in any of these cases.

System NaCl - Na₂SO₄ - KCl - K₂SO₄ - NaK₃(SO₄)₂ - Na₂SO₄.10H₂O - H₂0

As in the previous set of cases, the salts that precipitate in systems that have the Na⁺, K⁺, Cl^-, and ions in water depend on temperature. At the temperature of 373.15 K, there is experimental evidence for formation of the following five salts: NaCl, Na₂SO₄, KCl, K₂SO₄, and NaK₃(SO₄)₂. In addition to these five salts, a sixth salt, Na₂SO₄.10H₂O, can precipitate at 298.15 K. In the formulation of the cases discussed here, we assume the possibility of forming all these six salts and let the procedure find which of them, if any, precipitate(s) in a given case. In all cases, the specified initial amount of water is 55.50844 moles (1 kg), but the calculated final amount of water in the liquid phase is smaller than this amount when the hydrated salt, Na₂SO₄.10H₂O, precipitates. We compare our results with those of Nicolaisen et al. (1993) who represented the properties of these systems in diagrams whose axes are:

where m denotes an ion molality.

Temperature of 373.15 K

Table 5 summarizes the specifications and results of cases at the temperature of 373.15 K. In agreement with experimental data, the procedure does not predict the precipitation of Na₂SO₄.10H₂O in any of these cases. In Case 11, the procedure predicts the precipitation of Na₂SO₄ and NaK₃(SO₄)₂. During the course of calculations, pure NaK₃(SO₄)₂ is the first solid phase added to the system, later followed by the inclusion of the pure solid phase of Na₂SO₄. Using the ion mole numbers of the saturated aqueous solution reported in Table 5, it is possible to compute and in Eqs. 24 and 25.

In Case 12, as in Case 11, the procedure predicts the precipitation of two salts, with pure NaK₃(SO₄)₂ being the first solid phase added during the course of calculations. K₂SO₄ is the other salt predicted to precipitate with the specifications of Case 12. Case 13 shows a case in which NaCl and Na₂SO₄ precipitate, with the latter being the first pure solid phase added to the system by the algorithm. The procedure predicts the precipitation of three salts in the conditions of Case 14 and adds them to the system during computations in the following order: NaK₃(SO₄)₂, NaCl, and KCl. Using the ion mole numbers obtained in Case 14 (Table 5) to compute molalities and utilizing them in Eqs. 24 and 25, the values of and are, respectively, equal to 0.53 and 0.09.

Case 15 illustrates the action of the phase elimination procedure. The addition of pure solid phases during problem solution occurs in the following order: NaK₃(SO₄)₂, KCl, and K₂SO₄. It is interesting to observe that the first solid phase added to the system by the phase stability test (pure solid phase of NaK₃(SO₄)₂) is removed during the four phase splitting calculation. The final condition in Case 15 has two pure solid phases of and K₂SO₄ in equilibrium with the aqueous ionic phase.

There is excellent agreement with the results of Nicolaisen et al. (1993) in Cases 11, 13, and 14, but somewhat less satisfactory in Cases 12 and 15 because the region of K₂SO₄ precipitation is slightly bigger in our calculations, a discrepancy whose origin we were unable to identify. Phase splitting calculations typically required three iterations of the Murray method, with the exception of Case 21 that required from four to six iterations, depending on the number of phases present.

Temperature of 298.15 K

Table 6 shows specifications and results at 298.15 K. Cases 16 and 17 show situations with the simultaneous precipitation of two salts, while three salts precipitate in Cases 18 and 19.

Cases 19 and 20 have similar specifications but rather different results. In Case 19, three salts precipitate: Na₂SO₄, NaK₃(SO₄)₂, and Na₂SO₄.10H₂O. Only two of them precipitate in Case 20: NaK₃(SO₄)₂ and Na₂SO₄.10H₂O. In the solution of Case 20, the first solid phase added to the system was Na₂SO₄, followed by Na₂SO₄.10H₂O and later by NaK₃(SO₄)₂. During the phase splitting computations with one liquid and three solid phases, the procedure removed the solid phase of Na₂SO₄ before achieving the final configuration. Phase splitting calculations typically required three or four iterations of the Murray method, but convergence of the intermediate three phase configuration (liquid phase + Na₂SO₄ + Na₂SO₄.10H₂O) in Case 20 required eight iterations.

Using the results of Table 6 to compute molalities and then and in Eqs. 24 and 25, there is excellent agreement with the results of Nicolaisen et al. (1993) in all cases.

CONCLUSIONS

The new procedure proposed here to compute solid-liquid equilibrium in electrolyte systems that may form pure solid phases at given temperature, pressure, and global composition has been shown to be reliable and fast. The procedure combines three sub-procedures: a phase stability test, minimization of the Gibbs free energy with a stoichiometric formulation of the salt-forming reactions to compute phase splitting, and a phase elimination test. After the phase splitting calculation, in which a given system has a certain number of phases, the phase stability test establishes whether including an additional phase will reduce the Gibbs free energy further. The criteria used for phase stability may lead to the premature inclusion of phases that should be absent from the final solution. Two of the examples illustrated this possibility and showed that the phase elimination sub-procedure removes such phases. It is possible to use the procedure with several excess Gibbs free energy models for liquid phase behavior, but we chose a model previously used in the literature for comparison. Agreement with these literature results is, in general, excellent.

ACKNOWLEDGMENTS

The authors acknowledge the financial support of PETROBRAS and its affiliates in the early stages of this work.

NOMENCLATURE

Roman Letters

a_S activity of salt S in its pure solid phase activity of salt S in the pure trial solid phase a_ij activity of species i in phase j G Gibbs free energy m_I molality of ion I molar mass of water in g / mol n_c number of species n_i number of moles of species i n_i0 initial number of moles of species i n_ij number of moles of species i in phase j n_I number of ionic species n_p number of phases R Universal gas constant T absolute temperature x_i mole fraction of species i mole fraction of species i in the trial phase

Greek Letters

γ_i activity coefficient of species i based on the mole fraction scale and symmetric convention activity coefficient of ion I based on the mole fraction scale and unsymmetric convention activity coefficient of ion I based on the molal scale and unsymetric convention µ_i chemical potential of species i µ_ij chemical potential of species i in phase j chemical potential of species i in trial phase chemical potential of salt S in its pure trial solid phase chemical potential (molar Gibbs free energy) of formation of pure solid salt S at the system temperature chemical potential (molar Gibbs free energy) of formation of pure liquid water at the system temperature ν_i stoichiometric coefficient of species i (negative for reactants, positive for products) ν_ik stoichiometric coefficient of species i in reaction k θ_ij fraction of species i present in phase j Ω_G phase stability function ξ_k extend of reaction k

Subscripts

I ion k S salt k reference state at 298.15 K k

Superscript

∇ formation properties from literature(Nicolaisen et al., 1993)

(Submitted: November 2, 2009; Revised: January 13, 2010; Accepted: January 13, 2010)

• Allison, J. D., Brown, D. S., Novo-Gradac, K. J., MINTEQA2/PRODEFA2, a Geochemical Assessment Model for Environmental Systems: Version 3.0 User's Manual. EPA/600/3-91/021, USA Environmental Protection Agency (1991).
• Alvarez, V. H., Larico, R., Ianos, Y., Aznar, M., Parameter Estimation for VLE Calculation by Global Minimization: The Genetic Algorithm. Brazilian Journal of Chemical Engineering, 25, 409 (2008).
• Bale, C. W., Chartrand, P., Degterov, S.A., Eriksson, G., Hack, K., Ben Mahfoud, R., Melançon, J., Pelton, A. D., Petersen, S., FactSage Thermochemical Software and Databases. Calphad, 26, 189 (2002).
• Ballard, A. L., Sloan, E. D., The Next Generation of Hydrate Prediction. Part III. Gibbs Energy Minimization Formalism. Fluid Phase Equilib., 218, 15 (2004).
• Beahm, E. C., Toghiani, R. K., Weber, C. F., Comparative Calculations of Solubility Equilibria. Oak Ridge National Laboratory ORNL/TM-1999/127 (2000).
• Bonilla-Petriciolet, A., Bravo-Sanchez, U. I., Castillo-Borja, F., Frausto-Hernandez, S., Segovia-Hernandez, J. G., Gibbs Energy Minimization Using Simulated Annealing for Two-Phase Equilibrium Calculations in Reactive Systems. Chemical and Biochemical Engineering Quarterly, 22, 285 (2008).
• Bruin, S., Phase Equilibria for Food Product and Process Design. Fluid Phase Equilib., 160, 657 (1999).
• Cardoso, M. J. E. M., Calculation of Phase and Chemical Equilibrium in Electrolyte Systems. Ph.D. thesis, Technical University of Denmark (1989).
• Castier, M., Solution of the Isochoric-Isoenergetic Flash Problem by Direct Entropy Maximization. Fluid Phase Equilib. 276, 7 (2009).
• Castier, M., Rasmussen, P., Fredenslund, Aa., Calculation of Simultaneous Chemical and Phase Equilibria in Nonideal Systems. Chem. Eng. Sci., 44, 237 (1989).
• Cezac, P., A General and Applied Thermodynamic Model for the Characterisation of Saline Systems. Desalination, 175, 269 (2005).
• Chen, C. C., Toward Development of Activity Coefficient Models for Process and Product Design of Complex Chemical Systems. Fluid Phase Equilib., 241, 103 (2006).
• Cherif, M., Mgaidi, A., Ammar, M. N., Abderrabba, M., Furst, W., Modelling of the Equilibrium Properties of the System H₃PO₄H₂O: Representation of VLE and Liquid Phase Composition. Fluid Phase Equilib., 175, 197 (2000).
• Cherif, M., Mgaidi, A., Ammar, M. N., Abderrabba, M., Furst, W., Representation of VLE and Liquid Phase Composition with an Electrolyte Model: Application to H₃PO₄H2O and H₂SO₄H₂O. Fluid Phase Equilib., 194--197, 729 (2002).
• Chiavone, O., Rasmussen, P., Modeling of Salt Solubilities in Mixed Solvents. Braz. J. Chem. Eng., 17, 117 (2000).
• Corazza, M. L., Corazza, F. C., Cardozo Filho, L., Dariva, C., A Subdivision Algorithm for Phase Equilibrium Calculations at High Pressures. Braz. J. Chem. Eng., 24, 611 (2007).
• Ferrando, N., Lugo, R., Mougin, P., Coupling Activity Coefficient Models, Henry Constant Equations, and Equations of State to Calculate Vapor-Liquid and Solid-Liquid Equilibrium Data. Chem. Eng. Process., 45, 773 (2006).
• Floudas, C. A., Gounaris, C. E., A Review of Recent Advances in Global Optimization. J. of Global Optimization, 45, 3 (2009).
• Gil, V. M. S., Paiva, J. C. M., Computer Simulations of Salt Solubility. J. Chem. Educ., 83, 170 (2006).
• Huber, C., Filella, M., Town, R. M., Computer Modelling of Trace Metal Ion Speciation: Practical Implementation of a Linear Continuous Function for Complexation by Natural Organic Matter. Computers & Geosciences, 28, 587 (2002).
• Kulik, D., GEM-Selektor Code Package: History and Perspective in GEM Chemical Thermodynamic Modeling. GEMS-Workshop, November 19/20, Empa, Dübendorf, Switzerland (2007).
• Lin, Y., Thomsen, K., de Hemptinne, J.-C., Multicomponent Equations of State for Electrolytes. AIChE J., 53, 989 (2007).
• Lothenbach, B., Matschei, T., Glasser, F., Kulik, D., GEMS-Workshop, November 19/20, Empa, Dübendorf, Switzerland (2007).
• Messnaoui, B., Ouiazzane, S., Bouhaouss, A., Bounahmidi, T., A Modified Electrolyte-Uniquac Model for Computing the Activity Coefficient and Phase Diagrams of Electrolytes Systems. CALPHAD, 32, 566 (2008).
• Michelsen, M. L., The Isothermal Flash Problem. Part I. Phase Stability. Fluid Phase Equilib., 9, 1 (1982a).
• Michelsen, M. L., The Isothermal Flash Problem. Part II. Phase Split Calculation. Fluid Phase Equilib., 9, 21 (1982b).
• MINEQL++, http://www.mineql.com/, accessed on May 12 ( 2009).
  » link
• Montastruc, L., Azzaro-Pantel, C., Pibouleau, L., Domenech, S., Use of Genetic Algorithms and Gradient Based Optimization Techniques for Calcium Phosphate Precipitation. Chem. Eng. Process., 43, 1289 (2004).
• Murray, W., Numerical Methods for Unconstrained Optimization. Academic Press, London and New York (1972).
• Nagatani, G, Ferrari, J., Cardozo Filho, L., Rossi, C. C. R. S., Guirardello, R., Oliveira, J. V., Corazza, M. L., Phase Stability Analysis of Liquid-liquid Equilibrium with Stochastic Methods. Braz. J. Chem. Eng., 25, 571 (2008).
• Nicolaisen, H., Rasmussen, P., Sorensen, J.M., Correlation and Prediction of Mineral Solubilities in the Reciprocal Salt System (Na⁺, K⁺) (Cl^-, SO4^2-) - H2O at 0-100şC. Chem. Eng. Sci., 48, 3149 (1993).
• Pinho, S. P., Macedo, E. A., Representation of Salt Solubility in Mixed Solvents: A Comparison of Thermodynamic Models. Fluid Phase Equilib., 116, 209 (1996).
• Rahman, I., Das, A. K., Mankar, R. B., Kulkarni, B. D., Evaluation of Repulsive Particle Swarm Method for Phase Equilibrium and Phase Stability Problems. Fluid Phase Equilib., 282, 65 (2009).
• Simoni, L. D., Lin, Y. D., Brennecke J. F., Stadtherr, M. A., Reliable Computation of Binary Parameters in Activity Coefficient Models for Liquid-Liquid Equilibrium. Fluid Phase Equilib., 255, 138 (2007).
• Templeton, D. M., Ariese, F., Cornelis, R., Danielsson, L.-G., Muntau, H., van Leewen, H. P., Lobinski, R., Guidelines for Terms Related to Chemical Speciation and Fractionation of Elements. Definitions, Structural Aspects, and Methodological Approaches. (IUPAC Recommendations 2000), Pure Appl. Chem. 72, 1453 (2000).
• Thomsen, K., Aqueous Electrolytes: Model Parameters and Process Simulation. Ph.D. Thesis, Technical University of Denmark (1997).
• Thomsen, K., Rasmussen, P., Gani, R., Correlation and Prediction of Thermal Properties and Phase Behaviour for a Class of Aqueous Electrolyte Systems. Chem. Eng. Sci., 51, 3675 (1996).
• Toghiani, B., Lindner, J. S., Weber, C. F., Hunt, R. D., Modeling of Sulfate Double-salts in Nuclear Wastes. Oak Ridge National Laboratory ORNL/TM-2000/317 (2000).
• University of Colorado at Boulder Interactive Simulations, http://phet.colorado.edu/simulations/ sims.php?sim=Salts_and_Solubility, accessed on May 5 (2009).
  » link
• van der Lee, J., Themodynamic and Mathematical Concepts of CHESS, Technical Report Nr. LHM/RD/98/39. Ecole des Mines de Paris (1998).
• van der Lee, J., DeWindt, L., CHESS Tutorial and Cookbook, Updated for version 3.0, Report Nr. LHM/RD/02/13. Ecole Nationale Supérieure des Mines de Paris (2002).
• Villafáfila García, A., Thomsen., K., Stenby, E. H., Prediction of Mineral Scale Formation in Geothermal and Oilfield Operations Using The Extended UNIQUAC Model Part I. Sulfate Scaling Minerals. Geothermics, 34, 61 (2005).
• Villafáfila García, A., Thomsen., K., Stenby, E. H., Prediction of Mineral Scale Formation in Geothermal and Oilfield Operations Using The Extended UNIQUAC Model - Part II. Carbonate-Scaling Minerals. Geothermics, 35, 239 (2006).

*

To whom correspondence should be addressed

On leave from Escola de Química, Universidade Federal do Rio de Janeiro. This work was partially done at COPPE/UFRJ.

Publication Dates

• Publication in this collection
  14 Apr 2010
• Date of issue
  Mar 2010

History

• Received
  02 Nov 2009
• Reviewed
  13 Jan 2010
• Accepted
  13 Jan 2010