Gibbs Energy Minimization Model for Solvent Extraction with

Article pubs.acs.org/est

Cite This: Environ. Sci. Technol. 2019, 53, 7736−7745

Gibbs Energy Minimization Model for Solvent Extraction with Application to Rare-Earths Recovery Chukwunwike O. Iloeje,*,† Carlos. F. Jove ́ Coloń ,∥ Joe Cresko,§ and Diane J. Graziano‡ Energy Systems Division and ‡Decision and Infrastructure Sciences Division, Argonne National Laboratory, 9700 South Cass Avenue, Lemont, Illinois 60439-6903, United States ∥ Sandia National Laboratory, 1515 Eubank SE, Albuquerque, New Mexico 87123, United States § Advanced Manufacturing Office, US Department of Energy, 1000 Independence Avenue SW, Washington, DC 20585, United States Downloaded via UNIV OF SUSSEX on July 21, 2019 at 03:19:42 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

S Supporting Information *

ABSTRACT: The emergence of technologies in which rare-earth elements provide critical functionality has increased the demand for these materials, with important implications for supply security. Recycling provides an option for mitigating supply risk and for creating economic value from the resale of recovered materials. While solvent extraction is a proven technology for rare-earth recovery and separation, its application often requires extensive trial-and-error experimentation to estimate parameter values and determine experimental design configurations. We describe a modeling strategy based on Gibbs energy minimization that incorporates parameter estimation for required thermodynamic properties as well as process design for solvent extraction and illustrate its applicability to rare earths separation. Visualization analysis during parameter estimation revealed a linear relationship between the standard enthalpies of the extractant and respective organo-metal complexes, analogous to the additivity principle for predicting molar volumes of organic compounds. Establishing this relationship reduced the size of the parameter estimation problem and yielded good agreement between model predictions and reported equilibrium extraction data, validating the property estimates for the organic phase species. Design exploration and optimization results map the space of feasible solvent extraction column configurations and identify the set of optimal design parameter values that meet recovery and purity targets. disruptions,12,13 and value recovery will benefit the local manufacturing industry.12,14 Solvent-based extraction technologies are widely considered suitable for recovering high-purity fractions of REEs.15 These liquid−liquid extraction systems take advantage of the relative ability of solutes to distribute between immiscible aqueous and organic phases in equilibrium. The aqueous solvent promotes dissociation of ionic species, while the organic phase contains a complex-forming extractant that combines with the target ionic species to form organo-metal complexes, with high solubility in

1. INTRODUCTION Rare-earth elements (REEs) have unique magnetic, catalytic, and phosphorescent properties that significantly improve performance of a wide range of technologies.1 These technologies span aerospace, clean energy, telecommunications, electronics, transportation, defense, and other diverse applications. For these applications, certain REEs are considered to be critical materials 2−7 with important implications for supply security, economics, and disposal management. Recycling REEs provides an option for closing the material flow loop, diversifying supply sources, and creating economic value from the resale of recovered material. Closing the material loop will conserve the primary resource,8−11 diversifying sources will insulate the supply chain from © 2019 American Chemical Society

Received: Revised: Accepted: Published: 7736

March 21, 2019 May 31, 2019 June 3, 2019 June 3, 2019 DOI: 10.1021/acs.est.9b01718 Environ. Sci. Technol. 2019, 53, 7736−7745

Article

Environmental Science & Technology

in each phase represents the equilibrium criterion for a multiphase, multicomponent system under the constraints of fixed temperature and pressure.19,22 For such a system, the Gibbs minimization problem takes the following form:

the organic phaseand, conversely, negligible solubility in the aqueous phase.15,16 Depending on the type of extractant, this complex formation may involve solvation, anion exchange, or cation exchange. In solvation, the extractant displaces water molecules from the coordination sphere of a metal complex; in anion exchange, an anionic metal complex from the aqueous phase displaces an anion from the organic extractant; in cation exchange, the metal cation displaces a proton from the organic extractant.16 The 2-ethylhexyl hydrogen 2-ethylhexyl phosphonate (HEH-EHP, commonly known as PC88A) extractant considered in this study belongs to the cationic group and is considered effective for REE recovery.15 Eq 1 illustrates the overall cationic exchange reaction for hydrometallurgical extraction, where Mn+ represents the cation, (HA)2 indicates the PC88A extractant, and M(HA2)n represents the organometal complex. M

n+

eq

+ n(HA)2,org ↔ M(HA 2)n ,org + nH

+

min G =

∑ ∑ nipGi̅ p α

(2)

i

subject to Gi̅ = Gi̅ 0 + RT ln(ai)

(3)

ai = γimi

(4)

∑ nip = ni

(5)

i

∑ cj ,i = dj

(6)

i

(1)

nip ≥ 0

In practice, the physical and chemical similarities of REE mixtures make extraction difficult, needing multiple separation stages to meet design targets.12,17 The actual design of the separation process often requires extensive trial-and-error experimentation to characterize the extraction chemistry and to determine process configurations that meet separation and recovery targets. These factors motivate theoretical modeling strategies to predict extraction equilibria and identify optimal process configurations. Approaches for modeling the thermodynamics of liquid−liquid equilibrium in solvent extraction systems build on the theoretical approach proposed by Gibbs,18 but they differ on the choice of physical quantity and solution algorithms.19,20 Gibbs energy minimization (GEM) computes liquid−liquid equilibrium by determining the system’s minimum Gibbs free energy with respect to the amount of each species in each phase. Unlike equilibrium constant-based approaches, it does not require knowledge of the actual solubility and speciation reactions, but instead it uses thermodynamic information about the chemical species composing in each phase at equilibrium.21 Yet it has seen limited use in REE solvent extraction literature, because, in many cases, essential standard thermodynamic properties of key organic-phase species are not available in the open scientific literature. In this study, we implement an approachbased on GEMfor estimating missing standard thermodynamic properties of organic-phase species and predicting extraction equilibria in hydrometallurgical systems. This GEM approachbuilding on earlier work by Jové Colón et al.21 reduces the requirement for extensive regression of experimental data to retrieve parameters. We illustrate its applicability to the design, analysis, and optimization of REE separations, with emphasis on mixture compositions relevant to recycling nickel metal hydride (NiMH) batteries. The first part of this paper describes the Gibbs minimization and standard thermodynamic property estimation method, while the second part illustrates an application to rare-earths separation by solvent extraction, covering multiparameter sweep studies for process design exploration and optimization analysis.

(7)

∑ nizi = 0

(8)

i

i ∈ {1, 2, ..., N}; j ∈ {1, 2, ..., M}

where npi and G̅ pi represent the number of moles and partial molar Gibbs energy of species i in phase p, G̅ 0i is the molar Gibbs free energy of species i in its standard state, ai is the activity, a measure of the effective concentration of species i that accounts for any deviations from ideality, γi is the activity coefficient, mi represents molality, the number of moles of solute per kilogram of solvent, and zi is the charge of species i in the aqueous phase. Eq 4 defines activity as the product of the activity coefficient and molality. Eq 5 satisfies species balance, where ni is the total number of moles of species i in all phases. Eq 6 ensures elemental balance, where dj is the total number of moles of element j across all species, and cj, i represents the number of atoms of element j per molecule of species i. Eq 7 satisfies the non-negativity criterion, while eq 8 represents charge balance. Solving eqs 2−8 predicts the species composition and phase distribution at chemical equilibrium. To solve, we first define the appropriate thermodynamic models for the electrolyte and organic phases. 2.2. Electrolyte Phase Chemistry. The activity coefficient model for the electrolyte phase is based on the Pitzer formulation for the excess Gibbs free energy in ionic solutions.23 The Pitzer formulation represents the nonideal behavior by accounting for short- and long-range interactions between ions in solution. In this study, the Pitzer model is used to represent the HCl electrolyte interactions, with electrolyte parameters refitted from Holmes et al.24 The ionic standard state for the aqueous solute is a hypothetical 1 molal solution referenced to infinite dilution.25 Standard partial molal thermodynamic properties of the aqueous species in the electrolyte are computed using the Helgeson-Kirkham-Flowers equation of state (EoS).26−29 This EoS provides the standard partial molal property values as a function of pressure and temperature, accounting for solvation effects that arise from ion−solvent interactions. Its expressions for the electrostatic contribution derive from Born equations, while the structure effects are determined from semiempirical expressions of the form in eq 9, regressed from experimental data.30

2. METHODS 2.1. Gibbs Energy Minimization. The minimum of Gibbs free energy with respect to the molar quantities of each species

ΔEi̅ 0 = f (αi , P , T ) 7737

(9) DOI: 10.1021/acs.est.9b01718 Environ. Sci. Technol. 2019, 53, 7736−7745

Article

Environmental Science & Technology Here, E represents some partial molar quantity (volume, compressibility, or specific heat); P and T are pressure and temperature, respectively, and αi represents the regressed virial coefficients. These coefficients have been documented for several electrolyte and ionic species by Shock et al.30,31 We provided the coefficient values for the ionic species as inputs to the electrolyte phase model defined in Cantera, an open-source suite of object-oriented, constitutive modeling tools and property databases, used for equilibrium calculations.25 (See Supporting Information section S0-A for the input data file specifications. For further information on defining the phase model in Cantera, see Moffat and Jové Colón.25) We validate the electrolyte model by predicting water activity and osmotic pressure for HCl−H2O systems and comparing these predictions with data from Hamer and Wu.32 Figure 1 shows good agreement for the molality range of interest.

includes all other input parameters like molar volumes, temperature, pressure, as well as the vector of experimental feed composition data. Yexp is any set of experimentally determined equilibrium isotherms, and f represents the corresponding model predictions determined by solving the GEM problem (eqs 2−8). A typical equilibrium isotherm is the Distribution ratio, D, a ratio of the amount of solute in the organic phase to the amount left in the aqueous phase at equilibrium. For this problem, we formulate the residual function with log(D̃ ), because the log-function improved gradient scaling, with D̃ defined as nj ,org + ε D̃ = f (D , ε) = nj ,aq + ε (11) where nj,org and nj,aq are the organic and aqueous phase molar amounts of species j at equilibrium, and ε is a very small value used for numerical smoothening to avoid division-by-zero exceptions. Figure 2 illustrates the parameter-estimation

Figure 1. Electrolyte model validation. Predictions for activity and osmotic coefficients match data from Hamer and Wu. Deviations become significant as molality increases beyond 8 m, but the region of interest for this study is below 4 m, where the fit is very good.

2.3. Organic-Phase Chemistry. The organic phase comprises the organic extractant (HA)2, an inert diluent, and the organo-metal complexes with the stoichiometry M(HA2)3. The corresponding solute and solvent standard states are the pure solute and pure organic solvent, respectively. We model this phase as an ideal mixture, which removes the complexity of predicting activities for the organo-metal complexes (γi = 1 in eq 4). This assumes that, by representing the organic complex association equilibria, the mixing model indirectly captures any nonideal behavior.21,33−35 The ideal mixing model requires the standard Gibbs energy values (G̅ 0i ) as input. However, the standard Gibbs energies for the extractant and organo-metal species are not readily available in open literature, so we estimate them by regression to reported equilibrium data as described in the following section. 2.4. Standard Gibbs Parameter Estimation. The standard Gibbs parameter estimation for the organic-phase species involves finding the minimum on a multidimensional residual function surface defined by the least-squares difference of predicted and experimental isotherms. Eq 10 expresses the functional form of the least-squares objective.

Figure 2. Thermodynamic parameter estimation logic.

process. It starts with an initial guess for the standard molar Gibbs energy for the organic species, which we provide indirectly as the standard entropy (S0) and the standard enthalpy of formation (ΔH0f ). We assume a fixed dummy value for the absolute entropy data for the organic species (approximated from the correlation by Glasser et al.36). In this way, ΔH0 replaces ΔG0 as the explicit model variable. The extractant standard molar enthalpy (ΔHk) was initialized with a default value (for our case, we used that of tributyl phosphate21), while those for the organo-metal complexes were initialized as nΔHk, where n is the ratio of the extractantto-metal mole proportion in the complex from eq 1. The extractant molar volume was obtained from vendor data.37 Those for the organo-metal complex were initially set at 3 times this value but later were refined by refitting. The refitting results showed that equilibrium predictions were largely insensitive to the specified organo-complex molar volume. The prediction step follows, where, for each experimental data point, the model minimizes the Gibbs energy of the system to predict multiphase equilibrium composition. Then, it uses the result to minimize the least-squares function of eq 10. If the problem converges, the process exits with the current standard molar property values; otherwise, it proceeds to an update step and repeats the cycle. The modeling framework

Ndata NRE

min

2 ∑ ∑ (f (ΔGk0 , ΔGj0 , φi ,j ,k) − Y iexp ,j ) i

j

(10)

Here, i refers to experimental data points, j refers to the set of organo-metal species, and k refers to the organic extractant. ΔG0 = ΔH0 − TΔS0 is the vector of species Gibbs energy, expressed as a function of enthalpy (H) and entropy (S). φi,j,k 7738

DOI: 10.1021/acs.est.9b01718 Environ. Sci. Technol. 2019, 53, 7736−7745

Article

Environmental Science & Technology

Figure 3. Plots of least-squares residuals for the simple case of a Nd-HCl-PC88A system with two fit parameters: the standard molar enthalpy for the extractant (ΔHPC88A) and that for the complex species (ΔHNd). (a) Contour plot illustrating degree of fit between predicted and experimental39 values. (b) The locus of points with high degree of fit (function value where relative error