A Robust and Efficient Algorithm for Computing Reactive Equilibria in

Aug 27, 2014 - ABSTRACT: Prediction of liquid−liquid phase equilibria in reacting systems is ... to adopt a theoretical approach for computing equil...
0 downloads 0 Views 600KB Size
Download PDF
Article pubs.acs.org/IECR

A Robust and Efficient Algorithm for Computing Reactive Equilibria in Single and Multiphase Systems M. V. S. R. Ravi Kanth,† S. Pushpavanam,*,‡ Shankar Narasimhan,‡ and Murty B. Narasimha† †

Nuclear Fuel Complex, Department of Atomic Energy, Hyderabad, India 500062 Department of Chemical Engineering, Indian Institute of Technology, Madras, India 600036

‡

S Supporting Information *

ABSTRACT: Prediction of liquid−liquid phase equilibria in reacting systems is important in many applications such as reactive extraction. This problem poses several numerical challenges. These systems are governed by highly nonlinear algebraic equations and are plagued by the issue of nonconvergence of iterative algorithms. They exhibit strong sensitivity to the choice of the initial values or starting guesses for the variables. The nonconvergence stems primarily from the wide range of concentrations of various species at equilibrium. In this work, a new methodology for predicting thermodynamic equilibria of multiphase reacting systems is proposed. The mathematical formulation and solution are based on the use of the logarithms of the concentrations as the dependent variables. The proposed algorithm shows rapid convergence even when the initial guesses are far from equilibrium. The efficiency of the method is demonstrated by predicting the equilibrium concentrations of species in three systems of varying degrees of complexity.

1. INTRODUCTION Liquid−liquid extraction is a common process used for commercial-scale purification of important metals in the mining and nuclear industries. In this process, one of the phases is an aqueous phase containing the desired metal in solution, and the other is an organic/solvent phase. Agitated extractors such as mixer−settlers are commonly used mass-transfer equipment for extraction, because equilibrium is attained quickly in them. Reactive extraction is an important class of liquid−liquid extraction processes in which, in addition to extraction, several chemical reactions occur, forming ions and solute−solvent complexes. The solvent extraction of iron1 and the purifications of zirconium,2,3 uranium and plutonium,4 and neptunium and other actinides,5 among other processes, are carried out through reactive extraction. In systems containing metals that are hydrolyzable, such as zirconium, plutonium, and neptunium, the efficiency of the extraction process is dependent on the ionic strength, which, in turn, is a function of the acidity and the concentrations of various metal complexes.5 The performance of these processes can be predicted using models based on thermodynamic equilibrium. The key objective in developing such models is to compute the concentrations of various species, such as ions and metal complexes, in the different phases at thermodynamic equilibrium, in order to estimate the partition or distribution coefficient (D). The distribution coefficient is a function of several variables such as temperature, acidity, metal concentration, and solvent concentration. The experimentally determined value of D is applicable only to the specific experimental conditions used in its determination. It is thus necessary to conduct large numbers of experiments to establish the partitioning of metals over the entire operating range. This is a labor-intensive task. Empirical models for the evaluation of distribution coefficients are not recommended because they cannot capture the complex interactions of variables in the © 2014 American Chemical Society

system. Hence, there is a need to theoretically estimate this parameter from fundamental principles. It is therefore necessary to adopt a theoretical approach for computing equilibrium concentrations. Equilibrium concentrations in heterogeneous liquid−liquid systems can be estimated using either a nonstoichiometric or stoichiometric approach.6 The former technique determines the thermodynamic equilibrium by minimizing the Gibb’s free energy of the system. The latter uses reaction stoichiometry and solves mass-action equations satisfying mass balance and electroneutrality. This gives rise to nonlinear equations representing the phase and reaction equilibria with mass and charge balances as constraints. It can be shown that the two methods are equivalent. 6 In this work, we use the stoichiometric method, which requires only the equilibrium constants of various reactions as inputs. Computational packages are available for predicting the equilibrium concentrations of species in aqueous solutions based on both the nonstoichiometric and stoichiometric methods of formulation.6 These packages are mostly applicable for geochemical systems. A tailor-made approach is essential for predicting the concentrations of species such as zirconium, plutonium, and neptunium in the context of estimating liquid− liquid equilibria. Morel and Hering7 described a solution technique called the “Tableau” method based on the stoichiometric approach. In this method, the available species are classified as “components” and “secondary species”. Their concentrations are evaluated by sequentially solving two basic sets of equations: mass-balance and electroneutrality equations and reaction equilibria. The Received: Revised: Accepted: Published: 15278

July 2, 2014 August 22, 2014 August 27, 2014 August 27, 2014 dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

concentrations as δ. Among the NS species, we consider NC of the species as components, whose concentrations are denoted as δC. These NC components can be viewed as forming an independent basis,7 in that all other species can be represented as the products of reactions involving only the components and no component can be written as the product of a reaction involving only the other components.10 Determining the concentrations of NS species requires NS independent equations. If there are NR independent reactions (NR < NS) associated with the formation of each species, then the reaction equilibria provide NR equations. The number of components (NC) to be selected is NC = NS − NR. The mass conservation and electroneutrality equations of these NC components constitute the remaining set of independent equations. Morel and Hering7 explained that the selection of components should be such that the components can form all of the species present in the system. It is preferable to select the components as the dominant species in the given system, as this helps in obtaining good guesses for their concentrations. A generalized procedure for identifying the set of components among the available species is detailed in the Appendix. In these systems, various species participate in reversible chemical reactions. As a representative of all chemical/ionic reactions, considering the formation of species C as

concentrations of the component species are obtained by using the Newton−Raphson method on modified mass-balance and electroneutrality equations. The concentrations of secondary species are deduced from reaction equilibria. This approach suffers from several disadvantages. For example, it is prone to the problem of nonconvergence. In addition, it is sensitive to the choice of the initial or starting guesses for the variables8,9 and, consequently, requires a large number of iterations. This method requires a good candidate for the initial guess that should be close to the true solution for convergence. Ensuring convergence of the algorithm for a new system would hence be a challenge. Carrayrou et al.9 described an improved algorithm for solving the thermodynamic chemistry by overcoming this convergence problem. Their approach requires that additional constraints be specified to make the algorithm robust. Recently, Leal et al.6 proposed an alternative numerical technique for multiphase equilibrium calculations to improve the convergence of highly nonlinear multiphase systems. However, this technique is also sensitive to the initial guess. It is based on a projection procedure to avoid negative concentrations in the final converged solution. From the above summary, it is clear that the theoretical basis for calculating phase equilibria in the presence of reactions is well understood. However, there are several challenges in obtaining the numerical solutions to the governing equations.6,8,9 The nonlinear algebraic governing equations have to be solved iteratively. The convergence of the algorithms is strongly dependent on the initial guess. In this case, the challenge arises because the concentrations of various species span a wide range of values. This work focuses on the development of a numerical algorithm that is robust, namely, that will converge for a wide choice of initial guesses. In this article, we present a new methodology for solving thermodynamic equilibria of single-phase and multiphase reacting systems. It has the advantage of ensuring convergence with initial guesses that are far from equilibrium. The method does not impose any new constraints and eliminates the sensitivity of previous methods to the initial guess. We detail the mathematical formulation in section 2, before proceeding to the details of the numerical algorithm in section 3. We then apply the proposed methodology for predicting equilibria in section 4. Specifically, we discuss its applicability for modeling aqueous equilibria of the CaCO3−water system, as well as solvent extraction equilibria for the ferric chloride− tributyl phosphate (TBP) system. The robustness of the numerical technique is further demonstrated by applying the methodology to the highly complex aluminum−gallic acid system.8,9

(1)

a A + b B ↔ cC + d D

These reactions can involve species that are present in different phases and can be expressed as NS

∑ sjiδi ↔ 0,

for all reactions j = 1, 2, ..., NR

i=1

(2)

where sji is the stoichiometric coefficient of the ith species in the jth reaction and, in accordance with convention, it is positive for products and negative for reactants. At thermodynamic equilibrium, the rates of the forward and backward reactions are equal. The equilibrium constant (K) of the jth reaction is written as NS

Kj(T , P , δ) =

∏ ai s , ji

for all reactions j = 1, 2, ..., NR

i=1

(3)

where T and P are the temperature and pressure, respectively, of the system; Kj is the equilibrium constant of the jth reaction; and ai is the activity of the ith species, which can be expressed in terms of concentration (δi) as ai(T , P , δ) = δiγi(T , P , δ)

(4)

where γi is the activity coefficient of the ith species. The reactions for which equilibrium data are available might include reactant species that have not been selected as components. To estimate the concentrations of the unknown species, it is essential to express all reaction products as species that are formed stoichiometrically from species that have been selected as components. In other words, each reaction should contain only one species that has not been selected as a component, whereas all other species participating in the reaction should belong to the selected components. To accomplish this, reactions that are not expressed as speciesformation reactions involving only components are rewritten by taking linear combinations of the available reactions. This is possible because all reactions are independent. The equilibrium constants for the modified reactions (Kj*) are then evaluated.

2. MATHEMATICAL FORMULATION In this section, we describe the mathematical formulation required for computing the equilibrium concentrations of species in the context of liquid−liquid equilibria (LLE) with reactions. Adopting the definitions of the Tableau method,7 we consider a system containing NS species. A species can be a nonionized compound (water, solvent, salt, etc.), an ion, or a metal−organic complex. In a multiphase system, some or all of the species can be present in more than one phase. In this article, we restrict our considerations to multiphase reacting systems in which a species is present in only one phase. We denote the molar concentration of a species i (regardless of the phase in which it is present) as δi and the vector of species 15279

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

As a consequence, each species j is viewed as being formed from a reaction involving the set or subset of NC components. This allows one to write each species-formation reaction as

In this work, to have a generic algorithm, we used the charge balance to ensure electroneutrality and did not impose the proton balance.

NC

∑ νjiδi ↔ δj ,

3. NUMERICAL METHOD In this section, we describe the numerical method used to predict liquid−liquid equilibria, that is, to estimate unknown species concentrations at equilibrium. For this purpose, a total NS independent nonlinear algebraic equations formulated from (i) the set of reaction equilibria given by eq 6, (ii) the component balances given by eq 8, and (iii) the charge balance given by eq 9 have to be simultaneously solved. One of the challenges in this problem is that the concentrations of the various species span a large range of values, namely, from 101 to 10−14. This is one of the main causes for the numerical nonconvergence. This problem can be overcome by scaling the concentration variables using the logarithms of the concentrations. Using the logarithmic scale for concentration has other benefits, because it also restricts the solution to positive concentrations; avoids complexities due to the presence of extremely small concentrations of species (as low as 10−14); and avoids complexities due to large differences in the magnitudes of species concentrations, which is primarily responsible for nonconvergence. Taking the logarithms of both sides of eq 6 yields

for all noncomponent species j = 1, 2,

i=1

..., NR

(5)

where νji is the stoichiometric coefficient of the ith component in the jth species. By definition, component species (the basis set) cannot be represented as being formed from other components. Therefore, eq 5 is limited to the formation reactions of NR noncomponent species only. The equilibrium formation constant (K*j ) of the jth species from the resultant set of reaction equilibria expressed as species-formation reactions from components is aj K *j = , for all noncomponent species j = 1, 2, NC ∏i = 1 ai νji ..., NR

(6)

The superscript * is used to denote the fact that the species is now written as being formed from components. A generalized procedure for modifying the available reactions and evaluating Kj* from the available data on Kj is detailed in the Appendix. Typically, liquid−liquid extraction operations are carried out at constant temperature and pressure (mostly ambient). It is hence justifiable to assume that K*j and aj vary only with concentration depending on the ionic strength of the system. The reaction equilibria for all independent species-formation reactions constitute NR equations. NC additional equations in the form of mass and electroneutrality constraints are required, to determine the complete set of NS unknowns. Conservation of mass for the selected components is invoked to generate this remaining set of equations. The total mass of the ith component present in the system, Mi, is conserved, yielding the expression

NC

log10 K *j = log10 aj −

i=1

for all noncomponent species j = 1, 2, ..., NR

log10 K* = log10 a − A log10 a C

for all components i = 1, 2, ..., NC (7)

where Q is the volume of the phase in which species j is present and Mwi is the molecular weight of the ith component. Equation 7 can simplified by dividing throughout by Mwi as NS j=1

for all components i = 1, 2, ..., NC (8)

where mi is the total number of moles of the ith component. Imposing the constraint that the aqueous phase is an electrolytic solution that is electrically neutral gives

log10 K* = log10 δ + log10 γ − A log10 δC − A log10 γC (12)

NS

In eq 12, log10 δ is the vector of logarithmic species concentrations, that is, [log10 δ1, log10 δ2, ..., log10 δNS]T, and log10 γ is the vector of logarithmic activity coefficients of species; both are arranged in the order of noncomponent species followed by component species. log10 δC is the vector of logarithmic component concentrations and log10 γC is the vector of logarithmic component activity coefficients. Equation 12 is the transformed set of reaction equilibria and is expressed in terms of logarithmic species concentrations as variables.

∑ zjδj = 0 j=1

(11)

In eq 11, log10 K* is the column vector comprising the logarithms of the equilibrium constants. A is the stoichiometric coefficient matrix, which contains the stoichiometric coefficients νji. In the present system with NS species having NR formation reactions, the stoichiometric coefficient matrix has NS rows arranged in the order of noncomponent species followed by component species and NC columns. a is the vector of activities of species again in the order of noncomponent species followed by component species. Finally, aC is the vector of component activities. Equation 11 consists of a system of NS equations. The first NR equations correspond to the reaction equilibria of the noncomponent species and are considered for the subsequent equilibrium calculations. Equation 11 can be further expressed in terms of concentrations using eq 4 as

j=1

∑ νjiδjQ = mi ,

(10)

Equation 10 can be expressed in matrix form as

NS

∑ νjiδjQM wi = Mi ,

∑ νji log10 ai ,

(9)

where zj represents the charge on the jth species. One of the phases in the liquid−liquid system is typically aqueous. As a result, one must account for either of H+ or OH− as one of the components, along with H2O. It has been demonstrated that proton balance is analogous to electroneutrality.7,11 Therefore, either charge balance or proton balance (but not both) constitutes an independent equation. 15280

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

3.1. Modified Mass-Balance Equations. Unlike the Tableau method, which is a sequential method of solving equations, our proposed methodology aims to solve the reaction equilibria simultaneously with the mass and charge balances. This simultaneous approach results in the solution of the full equilibrium problem. This is an important difference between the proposed method and the Tableau-based sequential solution approach for predicting equilibria. The other important difference is that we use the concentrations of species in transformed form, namely, the logarithms of the concentrations, as the dependent variables. To accomplish this, it is necessary to modify the mass and charge-balance equations so that they have logarithmic concentrations as variables. This ensures that the same variables are present in all of the equations and enables their simultaneous solution. The mass balances, eq 8, are modified in terms of logarithmic concentrations as

The proposed approach is robust. It is very convenient to generate initial guesses with a minimal understanding of the physics of the problem. In our approach, guess values for only the component concentrations (δC0) are required. The component concentrations are relatively easier to guess because they are dominant, in accordance with the method used to select components. These initial guesses are used to generate initial logarithmic-scaled concentrations of all of the species (log10 δ0) using the reaction equilibrium expression10 log10 δ0 = log10 K* + A log10 δC0

The actual concentrations of species at equilibrium (δeq) are further derived from the solver output using δeq = 10 log10 δ

δQ 10 j

) = mi ,

for all components i = 1, 2,

j=1

..., NC

(19)

These equilibrium concentrations are useful for evaluating the distribution coefficient (D), which can be used for further verification of the thermodynamic model with experimental data. The key steps of the algorithm as shown in Figure 1 were implemented using MATLAB.

NS

∑ νji(10 log

(18)

(13)

Usually, unit volumes of the aqueous and organic phases (Q) are considered, so the above equation in matrix form simplifies to AT(10 log10 δ) = m

(14)

where 10 log10 δ = [10 log10 δ1, 10 log10 δ2 , ···, 10 log10 δ NS ]T

and m = [ m1 , m2 , ···, m NC ]T

is the vector of the total numbers of moles of the components. The charge balance, eq 9, is also transformed in terms of logarithmic concentrations as NS

δ 10 i

∑ zi(10 log

)=0

i=1

(15)

which can also be expressed in matrix form as z T(10 log10 δ) = 0

(16)

where z is the vector of species charges, [z1, z2, ···, zNS]T. The set of equations consisting of eqs 12, 14, and 16 is solved using a nonlinear solver such as the Newton−Raphson method. This solves the system of nonlinear equations by driving the function residuals to zero simultaneously. The function residuals are expressed as ⎡ log K* + A log δC + A log γ − log δ − log γ ⎤ 10 10 C 10 10 ⎥ ⎢ 10 ⎥ f = ⎢ AT(10 log10 δ) − m ⎢ ⎥ ⎢ z T(10 log10 δ) ⎥ ⎣ ⎦

Figure 1. Proposed methodology for determining equilibrium concentrations of species.

3.2. Convergence Criteria. The iterative procedure stops whenever the results of an iteration satisfy

(17)

|(log10 δ)n + 1 − (log10 δ)n |

The algorithm thus seeks to obtain log10 δ, such that f(log10 δ) =0 Equation 17 constitutes a set of (NS + 1) equations, but ignoring the proton balance from the component mass-balance equations ensures that the system of equations has only NS independent equations.

|(log10 δ)n + 1|

< εn (20)

and the function residual vector fulfils the condition:

|f n + 1| < εf 15281

(21)

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

where εn and εf are the tolerance values. A double-precision number format with 15 digits after the decimal point was used in all calculations. The εn and εf values were each set as 10−14.

the solution was assumed to be ideal, and the activity coefficients of the species were considered to be unity. The activity of water was also assumed to be unity, as it was present in a very large proportion. A demo version of the software MINEQL12 was used to estimate the equilibria of this system by assuming the concentration of water to be constant at its initial value. To simulate similar conditions and allow a valid comparison,12 the concentration variation of H2O was considered to be insignificant, using the proposed methodology as well. Thus, in this case, the concentration of H2O was not considered as a decision variable, and correspondingly, the mass balance for water was omitted from the set of equations. The equilibrium concentrations of species predicted using our proposed methodology under this assumption are consistent with the results of MINEQL,12 as shown in columns 7 and 8 of Table 2. In this work, the equilibrium concentrations of species were also predicted considering the variation in the concentration of H2O. The dominance of the species H2O can be seen from Table 2. It is interesting to note that, although the equilibrium concentration of water was almost the same as the initial value, the concentrations of most species differed by an order of magnitude, and in the case of carbonic acid, the concentration was found to be 2 orders of magnitude higher. This shows that even a small variation in the concentration of water can give rise to significant changes in the concentrations of other species. To assess the accuracy of the estimates, we checked the converged solution using the mass-balance expression. The total numbers of moles of Ca2+ and CO32− computed using eq 8 from the equilibrium concentrations were compared with initial values used for these species in the simulation (10−3 g· mol each). The deviations of the estimated total numbers of moles from what was fed was used as a measure of accuracy. In this context, the error in the predictions using the proposed methodology was obtained as 10−9%. We therefore conclude that the theoretical equilibrium concentrations computed using our methodology were accurate. The model was further tested for efficiency in convergence with different choices of initial guesses. As shown in Figure 2, convergence was obtained, with the defined tolerance (10−14), for various combinations of initial guesses of Ca2+ and CO32− concentrations, keeping the concentration of H2O constant as 55.5556 M and the guess concentration of H+ as 10−10 M. The rate of convergence was established from the reduced number of iterations for convergence even with initial guesses that were far from equilibrium. 4.2. Solvent Extraction Equilibria of FeCl3 with TBP.1 To test the proposed methodology for estimating equilibria in a multiphase system, we selected the system of solvent extraction of ferric chloride (FeCl3) present in aqueous HCl solution being extracted by tributyl phosphate (TBP) solvent. This choice was motivated by the fact that this is a well-understood system.1 The main intention was to model the reactive extraction system (FeCl3 with TBP) without a constraint on the initial guesses for predicting liquid−liquid equilibria. Here, we used the reaction equilibria and thermodynamic values from Lee et al.,1 as indicated in Table S1 of the Supporting Information. At equilibrium, the aqueous phase was assumed to contain the following 13 species: FeCl2+, FeCl2+, FeCl30, FeCl4−, FeOH2+, Fe(OH)2+, Fe(OH)30, Fe2(OH)24+, OH−, Fe3+, Cl−, H+, and H2O. The organic phase was assumed to contain the following 3 species: FeCl3·HCl·2TBP, TBP, and toluene. It

4. CASE STUDIES In this section, we discuss the results obtained by applying the proposed solution methodology to three different systems. We first demonstrate the method to estimate the aqueous equilibria in a classical electrolytic system, namely, CaCO3−water. Thereafter, we show how this methodology is useful for estimating equilibria in a multiphase liquid−liquid system: the ferric chloride solvent extraction system with tributyl phosphate (TBP). We finally demonstrate the robustness of the methodology by applying it to the aluminum−gallic acid system. The latter is known to be prone to nonconvergence problems and is sensitive to the choice of initial guess.8,9 The efficiency of our proposed method is thus illustrated on this system. 4.1. CaCO3−Water System. To detail the proposed methodology, we first chose a single-phase aqueous electrolytic system resulting from the dissolution of CaCO3 in water. The main intention was to model the CaCO3−water system and predict its behavior without having to make an accurate initial guess for the equilibrium concentrations. Here, we use the reaction equilibria and thermodynamic values from the demo version of MINEQL+ 4.6,12 as indicated in Table 1. Table 1. Reaction Equilibria of the CaCO3−Water System at 25 °C12 a log10 K

reaction

H+ + OH− ↔ H 2O Ca

2+

+

+ CO3

H + CO3

2−

+

2H + CO3 Ca

2+

Ca

2+

2−

+

↔ CaCO3

3.21

HCO3−

10.3

↔

2−

14

16.7

↔ H 2CO3

+ H + CO3

2−

↔ CaHCO3 +

+ H 2O ↔ CaOH + H

+

+

11.6 −12.7

a

The following species (NS = 10) are considered to be present at equilibrium: CaCO3, HCO3−, H2CO3, CaHCO3+, CaOH+, OH−, Ca2+, CO32−, H+, and H2O.

This system involves 10 species (NS) with 6 independent reactions (NR), so the number of components to be chosen is 4 (NC = NS − NR). We selected Ca2+, CO32−, H+, and H2O as the components. It can be proved that the final solution for equilibrium concentrations is independent of the choice of components.7 The set of reactions was initially expressed as speciesformation reactions from the selected components using the generalized procedure described in the Appendix. The procedure detailed in the Appendix was used to obtain the stoichiometric matrix A and the equilibrium constant vector (log10 K*). The elements of A and the equilibrium constant values for the formation reactions of noncomponent species are provided in Table 2. The proposed procedure was used to estimate the equilibrium concentrations of species produced upon dissolving 0.001 g-mol of CaCO3 in 1 L of water. The resulting dilute electrolytic solution has a very low ionic strength. Therefore, 15282

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

Table 2. Stoichiometric Coefficient Matrix, Equilibrium Constants, and Equilibria for the CaCO3−H2O System equilibrium concentration estimated by Ca2+

a

CO32−

H+

H2O

log10 K*

CaCO3 HCO3− H2CO3 CaHCO3+ CaOH+ OH−

1 0 0 1 1 0

1 1 1 1 0 0

0 1 2 1 −1 −1

0 0 0 0 1 1

3.21 10.3 16.7 11.6 −12.7 −14

Ca2+ CO32− H+ H2O total (M)

1 0 0 0 10−3

0 1 0 0 10−3

0 0 1 0

0 0 0 1 55.5556

0 0 0 0

MINEQL12

proposed methodology excluding H2O

Species 10−4 10−4 10−8 10−6 10−6 10−4 Components 6.30 × 10−4 3.56 × 10−4 3.64 × 10−11 a

3.64 2.77 2.27 3.24 3.47 2.77

× × × × × ×

3.63 2.69 2.51 3.39 3.39 2.69

× × × × × ×

proposed methodology including H2O

10−4 10−4 10−8 10−6 10−6 10−4

8.79 8.33 1.39 1.47 1.47 8.36

6.30 × 10−4 3.64 × 10−4 3.71 × 10−11 a

× × × × × ×

10−5 10−4 10−6 10−5 10−5 10−4

8.83 × 10−4 6.28 × 10−5 6.65 × 10−10 55.5547

Considered to be invariant.

new set of reaction equilibria, the stoichiometric coefficient matrix (A) and equilibrium constant vector (log10 K*) were then obtained as detailed in Table S2 of the Supporting Information. Activity coefficients of the solutes in the aqueous phase were calculated using the Bromley equation1,13 log10 γM = −

0.5108(z M)2 I 0.5 1 + I 0.5

+ FM

(22)

where ⎡ ⎤ ⎢ (0.06 + 0.6B )|z z | ⎥ (|z | + |z |)2 X MX M X ⎥ M FM = ∑ ⎢ + B [X] MX 2 4 1.5 ⎢ ⎥ X 1 + I ⎢⎣ ⎥⎦ |zMzX|

(

)

(23)

and NS

I = 0.5 ∑ zi 2δi

Figure 2. Convergence map for various initial guesses for Ca2+ and CO32− with a constant H2O concentration and fixed guess H+ concentration. Convergence was obtained for all initial guesses: logarithms of the concentrations of Ca2+ and CO32−. The dark circle indicates the final equilibrium concentrations (6.30 × 10−4 M Ca2+ and 3.64 × 10−4 M CO32−). The number of iterations required for convergence is shown as a color bar.

(24)

i=1

γM is the activity coefficient of cation M, zM is the ionic charge of the corresponding ion, I is the ionic strength of the medium, and BMX is the interaction parameter between cation M and anion X. [X] represents the concentration of X. The activity coefficients of neutral species and species present in the organic phase were assumed to be unity. This assumption is justified because such species have insignificant ionic interactions and no charge. All parameters, including the Bromley interaction parameters, were obtained from Lee et al.,1 and the activity of water (aH2O) was calculated using the Pitzer equation14 for a strong electrolyte/ionic medium (NX)

should be noted that HCl was not considered as a separate species (even though it was one of the reactants) because it was assumed to be completely ionized. There are thus NS = 16 species, with each being present in either the aqueous phase or the organic phase, but not both. Because there are 10 independent species-formation reactions (NR), the number of components to be chosen is 6 (NC = NS − NR). These were selected as Fe3+, Cl−, H+, H2O, TBP, and toluene. Because none of the species are distributed between the two phases, there is no need to include any phase equilibrium equation, and only reaction equilibrium equations for noncomponent species and mass balances for the component species need to be solved simultaneously. The set of reactions was initially expressed as speciesformation reactions from the selected components using the generalized procedure described in the Appendix. Using this

log10 a H2O = −

2ϕ[NX] 55.51 ln 10

(25)

where Φ is the osmotic coefficient in a strong ionic medium and [NX] represents the concentration of the ionic medium NX. The resulting system of equations was solved using the proposed methodology, considering unit volumes of the aqueous and organic phases. Convergence was found to be rapid and efficient. The effects of variations in the HCl, FeCl3, 14

15283

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

and TBP concentrations on the extraction of FeCl3 were investigated by individually varying the initial concentrations of HCl from 0.5 to 3.5 kmol/m3, of FeCl3 from 0.1 to 1.5 kmol/ m3, and of TBP from 1 to 3 kmol/m3. Table S3 (Supporting Information) indicates that the distribution coefficient of iron (DFe) increased with increasing concentrations of HCl, FeCl3, and TBP. The predictions obtained using the present methodology are comparable to the published data1 and exhibit good agreement with the experimental data, as shown in Figure 3. In our system, iron does not exist in both phases in the same

Figure 4. Convergence map for various initial guesses of Fe3+ and Cl− concentrations. The dark circle is the equilibrium value (0.005 kmol/ m3 Fe3+ and 2.25 kmol/m3 Cl−) for the initial concentrations of 1 kmol/m3 FeCl3, 2.5 kmol/m3 HCl, and 3 kmol/m3 TBP. Convergence was obtained for all initial guesses, even those far from equilibrium. The number of iterations required for reaching the final equilibrium is shown as a color bar.

tions of this complex aluminum−gallic acid system using a wide range of initial guesses. We adopt the stoichiometric coefficient matrix, definitions of species and components, and other required data from Carrayrou et al.9 In this system with NS = 18 species and NR = 14 independent reactions, the number of components chosen is 4 (Al3+, H3L, H+, and H2O). The species and components of the system are listed in Table S4 of the Supporting Information. The activity coefficients of neutral species were assumed to be unity, and those of ionic species were calculated using the Davies equation.9 For comparison with the previously published results of Carrayrou et al.,9 we assumed that the precipitation of gibbsite does not occur and that water is present in excess and has a constant composition. Thus, the concentration of water was not a decision variable, and the material balance for water was not included as an equation. Carrayrou et al.9 also reported the equilibrium compositions of the species at a fixed equilibrium pH of 5.8. To simulate similar conditions, the H+ ion concentration at equilibrium was fixed corresponding to this pH value. Because the proton composition was known, the proton balance/ electroneutrality condition was omitted as a redundant equation. Consequently, we solved a system of 16 algebraic equations. The proposed numerical method was employed to solve the resultant system of equations. Convergence was again found to be rapid and efficient. The model output was comparable to the published data,9 as shown in Table S4 of the Supporting Information. The accuracy of the estimates was established through a check using the mass-balance expression in eq 8. The total numbers of moles of Al3+ and H3L used for simulation were 10−3 g·mol each. The tolerance or deviation of the estimated total numbers of moles from eq 8 compared to the actual numbers of moles was used as a measure of accuracy. The equilibria reported by Carrayrou et al.9 were found to be

Figure 3. Comparison of distribution coefficients of iron (DFe) calculated using the proposed methodology with experimental and published data.1

form, and consequently, the definition of the distribution coefficient of iron (DFe) is considered as (iron concentration in the organic phase) (iron concentration in the aqueous phase) [FeCl3·HCl · 2TBP] = [iron present in the aqueous phase in all forms]

DFe =

(26)

The proposed methodology was also tested with random values of initial guesses far from equilibrium (i.e., far from the solution). The algorithm was able to converge to the solution rapidly, within the defined tolerance, even with random values of initial guesses as indicated in Figure 4. The accuracy of the estimates was further established through a mass-balance check using eq 8. The total numbers of moles of FeCl3, HCl, and TBP used for the simulation were 1, 2.5, and 3 kmol, respectively. A deviation of 10−14 kmol was obtained between the estimated total numbers of moles from eq 8 and the actual total numbers of moles, which indicates that the proposed methodology was effective in obtaining accurate estimates. 4.3. Aluminum−Gallic Acid System.9 Brassard and Bodurtha8 and Carrayrou et al.9 both reported that the aluminum (Al)−gallic acid (H3L) system is highly nonlinear and has issues of nonconvergence. Specifically, this complex system is highly sensitive to initial guesses. In this section, we demonstrate the robustness of the proposed methodology by solving the equilibrium composi15284

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

accurate to within 10−4 g·mol as compared to 10−14 g·mol using the proposed methodology in this work. These results clearly indicate the capability of the proposed methodology in obtaining accurate estimates. Brassard and Bodurtha8 presented the convergence map illustrating the initial-guess sensitivity of this system. They reported that, for many starting points of logarithmic concentrations log [H3L] and log [Al3+] from −12 to −3, the Newton−Raphson solver becomes caught in an infinite loop. As a result, nonconvergence persists with several combinations of the initial guesses, and successful convergence was obtained by the authors only for (i) a narrow band around log [Al3+] = −4 with any value of log [H3L] and (ii) values of −7 < log [H3L] < −3 and log [Al3+] > −11. It is also noticed that the results of the convergence map for the Al−H3L system reported by Carrayrou et al. (Figure 1)9 were the converse of what was actually presented by Brassard and Bodurtha.8 We have tested the initial-guess sensitivity of the proposed methodology on this complex system. Using the same simulation conditions as described by Brassard and Bodurtha8 and Carrayrou et al.,9 convergence was obtained for all initial guesses, as indicated in Figure 5. It was also found that, with the

deduced from phase and reaction equilibria, satisfying the mass and charge balances. Use of logarithmic concentrations as decision variables in the governing equations (modified mass balances, charge balances, and reaction equilibria) is the unique feature of this proposed methodology. The efficiency of the proposed methodology in predicting equilibrium concentrations of species in highly complex liquid− liquid systems was demonstrated successfully. The implementation of this methodology increased the robustness of the algorithm. This algorithm converges for a wide choice of initial guesses, even those that are far from equilibrium. In addition to elimination of the initial-guess dependency, this methodology avoids converging to negative concentrations as solutions, without the need for any new constraints. Another important advantage is the rapidity or relatively low computation time with quick convergence, as the number of iterations required for convergence is minimal. Moreover, the numerical accuracy of the solutions obtained using the methodology proposed in this work is far superior to that obtained using methods presented by earlier researchers.

■

APPENDIX In this appendix a formal framework for identifying a set of components from a given set of species is presented. We also discuss a method by which a noncomponent species can be written as being formed from the set of identified components. A.1. Generalized Procedure for Identifying the Set of Components

Let E be the atomic matrix (NE × NS) for a system with NS species formed from NE elements. Let RE be the rank of E. Then, the maximum number of independent reactions (NR) that can be written for this system is NS − RE.15 If S = {sji} is the stoichiometric matrix (NR × NS) for the given set of reactions (for which equilibrium data are known), then S must be in the null space of E {i.e., ES = 0, which also implies that the maximum number of independent reactions is dim[N(E)] = NS − RE}. One can choose a set of independent components such that their columns in E are linearly independent (the maximum number of independent columns one can choose is RE). A.2. Generalized Method for Obtaining a Set of Reactions Such That Each Reaction Contains Only One Noncomponent Species and Deducing the Respective Equilibrium Constants

Figure 5. Convergence map obtained using the proposed methodology with various initial guesses. Convergence was obtained for all initial guesses: logarithms of the concentrations of Al3+ and H3L. The black circle indicates the final equilibrium (2.59 × 10−5 M Al3+ and 2.06 × 10−7 M H3L). The number of iterations required for convergence is shown as a color bar.

Taking the logarithm of both sides of eq 3 yields NS

log10 Kj =

∑ sji log10 ai ,

for all reactions j = 1, 2,

i=1

..., NR

(27)

proposed methodology, the rate of convergence was high. The number of iterations required for convergence within the defined tolerance level was found to be low.

Equation 27 can be expressed in matrix form as

5. SUMMARY AND CONCLUSIONS In this work, a robust and efficient algorithm is proposed for predicting equilibrium concentrations in liquid−liquid systems. These systems are highly nonlinear and are plagued by the issue of nonconvergence when the classical Tableau method based on the sequential method of solution is used. The methodology proposed in this work is based on simultaneously solving the nonlinear system of equations

The columns of stoichiometric matrix S are to be permuted such that the first NC columns correspond to the noncomponent species and the rest of the columns correspond to the component species. S = [SNC SC].

log10 K = S log10 a

log10 K = [ S NC SC ] log10 a

(28)

(29)

Note that SNC is a square invertible matrix (if a maximum set of independent reactions has been specified) 15285

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286

Industrial & Engineering Chemistry Research

Article

S NC−1 log10 K = [ I S NC−1SC ] log10 a

(30)

log10 K* = S* log10 a

(31)

methyl phosphine oxide and tributyl phosphate. Solvent Extr. Ion Exch. 1998, 16, 1047−1066. (4) Benedict, M.; Pigford, T. H; Levi, H. W. Nuclear Chemical Engineering, 2nd ed.; McGraw-Hill: New York, 1981. (5) Sarsfield, M. J.; Sims, H. E.; Taylor, R. J. Modelling Np(IV) Solvent Extraction between 30% Tri-Butyl Phosphate and Nitric Acid in the Presence of Simple Hydroxamic Acids. Solvent Extr. Ion Exch. 2011, 29, 49−71. (6) Leal, A. M. M.; Blunt, M. J.; LaForce, T. C. A robust and efficient numerical method for multiphase equilibrium calculations: Application to CO2−brine−rock systems at high temperatures, pressures and salinities. Adv. Water Resour. 2013, 62, 409−430. (7) Morel, F. M. M.; Hering, J. G. Principles and Applications of Aquatic Chemistry; Wiley-Interscience: New York, 1993. (8) Brassard, P.; Bodurtha, P. A feasible set for chemical speciation problems. Comput. Geosci. 2000, 26, 277−291. (9) Carrayrou, J.; Mośe, R.; Behra, P. New efficient algorithm for solving thermodynamic chemistry. AIChE J. 2002, 48, 894−904. (10) Batley, G. E. Trace Element Speciation: Analytical Methods and Problems; CRC Press: Boca Raton, FL, 1989. (11) Brezonik, P. L.; Arnold, W. A. Water Chemistry: An Introduction to the Chemistry of Natural and Engineered Aquatic Systems; Oxford University Press: New York, 2011. (12) MINEQL+ 4.6 Demo Version; Environmental Research Software: Hallowell, ME, 2007; available at http://www.mineql.com (accessed Aug 22, 2014). (13) Bromley, L. A. Thermodynamic properties of strong electrolytes in aqueous solutions. AIChE J. 1973, 19, 313−320. (14) Zemaitis, J. F., Jr.; Clark, D. M.; Rafal. M.; Scrivner, N. C. Handbook of Aqueous Electrolyte Thermodynamics: Theory & Application; Wiley: New York, 1986. (15) Reklaitis, G. V. Introduction to Material and Energy Balances; John Wiley & Sons: New York, 1983.

where S* = S NC−1S = [ I S NC−1SC ] = [ I S*C ]

(32)

S*C = S NC−1SC

(33)

S* is the desired set of reactions because each noncomponent species appears in only one of the modified reactions. K* represents the equilibrium constants of the formation reactions for the species from the selected components. Comparing eqs 30 and 31, the logarithm of the equilibrium constants for the new reactions can be computed as log10 K* = S NC−1 log10 K

(34)

NR

log10 K*j =

∑ S NC,ji−1 log10 K i, i=1

for all reactions j = 1, 2, ..., NR

SNC,ji−1

(35)

S−1NC.

where represents the elements of The stoichiometric matrix denoted as A in eq 11 is given by

■

⎡−S −1S ⎤ A = ⎢ NC C ⎥ ⎣I ⎦

(36)

ASSOCIATED CONTENT

S Supporting Information *

Reaction equilibrium details for the FeCl3−TBP system (Table S1); stoichiometric coefficient matrix and equilibrium constants for the FeCl3−TBP system (Table S2); experimental and calculated values for the distribution coefficient of iron (Table S3); and stoichiometric coefficient matrix, equilibrium constants, and results of equilibria estimated for the Al−H3L system at a fixed pH of 5.8 (Table S4). This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +91-44-22574161. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS M. V. S. R. Ravi Kanth and M. B. Narasimha wish to express their sincere gratitude to Shri N. Saibaba, Chief Executive, NFC, for providing the opportunity, support and permission for carrying out this work. Thanks are also due to Shri R. N. Jayaraj, Former Chief Executive, NFC, Smt Meena Ravindran, GM, Shri H. R. Ravindra, DGM, and Shri Johnson D’ Souza, SM, for their constant encouragement and support.

■

REFERENCES

(1) Lee, M. S.; Lee, G. S.; Sohn, K. Y. Solvent Extraction Equilibria of FeCl3 with TBP. Mater. Trans. 2004, 45, 1859−1863. (2) Alcock, K.; Bedford, F. C.; Hardwick, W. H.; McKay, H. A. C. Tri-n-butyl phosphate as an extracting solvent for inorganic nitrates I: Zirconium nitrate. J. Inorg. Nucl. Chem. 1957, 4, 100−105. (3) Brewer, K. N.; Herbst, R. S.; Todd, T. A.; Christian, J. D. Zirconium extraction into octyl (phenyl)-N,N-diisobutylcarbamoyl15286

dx.doi.org/10.1021/ie502639a | Ind. Eng. Chem. Res. 2014, 53, 15278−15286