FUNDAMENTAL STUDIES AND ENGINEERING MODELING OF HYDROGEN BONDING By Aseel Mohamed Ahmed Bala Ahmed A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Chemical Engineering - Doctor of Philosophy 2018 ABSTRACT FUNDAMENTAL STUDIES AND ENGINEERING MODELING OF HYDROGEN BONDING By Aseel Mohamed Ahmed Bala Ahmed This project aims to enhance the engineering modeling of hydrogen bonding, or association, by blending ab initio quantum calc ulations, fundamental molecular level findings from experimental techniques, and thermodynamic models . B ecause of the ubiquity of hydrogen bonding, applications for an improved association model are extensive, ranging from drug design to plastics manufactu ring. Therefore, a substantial amount of work has been aimed at improving traditional thermodynamic tools, which often fail to capture the behavior of associating systems accurately. To guide models, spectroscopic techniques have been leveraged to gain ins ight into the interactions between molecules in the liquid phase, but interpretation is difficult. Moreover, with the advancement of computational chemistry technology, molecular dynamics (MD) and quantum mechanical (QM) calculations have also been utilize d to understand the characteristics of hydrogen bonded clusters. However, few studies have combined all 3 techniques (the thermodynamic model, spectroscopy and ab initio calculations) in a rigorous way. To this end, an activity coefficient model for associ limitations are explored with parameters from literature. Furthermore, a sequential MD and QM protocol is designed which facilitates the interpretation of the hydroxyl vibrati on in infrared spectroscopy and a method is developed to quantify the entire band. Finally, the methods are used to calculate the value of the association constant for an alcoh ol + alkane system . iii This dissertation is dedicated to the memory of my mother , Amal Ziada, who gave everything and then some . iv ACKNOWLEDGEMENTS The work documented in this dissertation would not have been possible without the contributions and support of the following individuals. First and foremost, I owe my deepest gratitude to m y advisor, Dr. Carl Lira, whose passions for thermodynamics and teaching have been a tremendous inspiration to me. His guidance, encouragement and patience were absolute, and I am thankful to have benefitted from his mentorship throughout my time as a grad uate student. I would also like their valuable insight, assistance and contributions. I acknowledge with gratitude the National Science Foundation for funding the proj ect as well as our collaborators in industry, Drs. Tim Frank, Paul Mathias, Navin Patel, Eric Cheluget, Dung Vu, You Peng, and Suphat Watanasiri, for their input and advice. I am also grateful to have worked with my colleagues and friends, William G. Killi an, Jackson Storer, Renming Liu and Dr. Lars Peereboom, who have assisted with data collection and contributed greatly to this work. In my time at Michigan State University, I have been fortunate enough to meet many wonderful people who became my support system, sounding boards and friends. I would like to thank Amrita Oak, Souful Bhatia, Annalisa Grunwald, Iman Nezam, Neda Rafat, GracieLou Klinger and Aritra Sarkar for the many laughs, sound advice and great company. Finally, I would like to express my s incerest gratitude for the life - long love and encouragement of my parents, Amal Ziada and Mohamed Ahmed Bala, my siblings, Asem, Leena, and Sarya, and my aunt, Iman Ziada. To this list, I add the friends that have played just as great a role in this journe y: Sarah Gard and Fatima El - Hussain. I am thankful for their kindness, compassion and unwavering support throughout the many years of our friendship. v TABLE OF CONTENTS LIST OF TABLES ................................ ................................ ................................ ....................... viii LIST OF FIGURES ................................ ................................ ................................ ........................ x KEY TO SYMBOLS AND ABBREVIATIONS ................................ ................................ ......... xv Chapter 1. In troduction to Thermodynamic Modeling of Associating Systems ............................ 1 1.1. Background ................................ ................................ ................................ .......................... 3 1.1.1. Equations of state (EOS) ................................ ................................ ............................... 4 1.1.2. Activity coefficient models ................................ ................................ ........................... 6 1.1.3. Association modeling ................................ ................................ ................................ .... 8 1.2. Resear ch objectives ................................ ................................ ................................ .............. 9 1.3. Dissertation outline ................................ ................................ ................................ ............ 10 Chapter 2. Relation of Concentration - based Equilibrium Constants to Activity - based Equilibrium Constants in the Liquid Phase ................................ ................................ ................................ ....... 12 2.1. Introduction ................................ ................................ ................................ ........................ 12 2.2. Gibbs energy and the activity - based equilibrium constant ................................ ................ 13 2.3. Gas phase activity - based equilibrium constant ................................ ................................ .. 14 2.4. Liquid phase activity - based equilibrium constant ................................ ............................. 15 APPENDICES ................................ ................................ ................................ .......................... 19 Appendix A: Development of the activity - based equilibrium constant for athermal solutions ................................ ................................ ................................ ................................ ............... 20 Appendix B: Effect of mass density on volume of mixing ................................ ................... 21 Chapter 3. Relation of Wertheim Association Constants to Concentration - based Equilibrium Constants for Mi xtures with Chain - forming Components ................................ ............................ 22 3.1. Introduction ................................ ................................ ................................ ........................ 22 3.1.1. Motivation ................................ ................................ ................................ ................... 22 3.1.2. Historical perspective ................................ ................................ ................................ .. 23 3.2. Theory and analysis ................................ ................................ ................................ ........... 25 3.2.1. Association term ................................ ................................ ....... 26 - based equilibrium constant ................................ ...................... 29 3.2.3. Relationship between , ................................ ................................ ............ 30 3.2.4. Association term ................................ ............. 32 ................................ ....... 34 3.2.6. Calculating association parameters ................................ ................................ ............. 36 ................................ ...... 38 3.3. Conclusions ................................ ................................ ................................ ........................ 40 APPENDICES ................................ ................................ ................................ .......................... 41 Appendix C: Derivation of the Wertheim association contribution to the activity coefficient ................................ ................................ ................................ ................................ ............... 42 vi Appendix D: Relationship between the fugacity coefficient and true compressibility factor ................................ ................................ ................................ ................................ ............... 44 Appe concentration of species ................................ ................................ ................................ ........ 45 Chapter 4. Relation of Wertheim and Chemical Theories in Systems with Association and S olvation ................................ ................................ ................................ ................................ ....... 47 4.1. Introduction ................................ ................................ ................................ ........................ 47 4.2. Theory and analysis ................................ ................................ ................................ ........... 50 4.2 ................................ ................... 50 4.2.2. Implementation of chemical theory ................................ ................................ ............ 52 4.2.3. Numerical equivalence ................................ ... 56 4.2.4. Activity coefficients ................................ ................................ ................................ .... 64 4.2.5. Perspective on further generalization ................................ ................................ .......... 66 4.3. Conclusions ................................ ................................ ................................ ........................ 67 APPENDICES ................................ ................................ ................................ .......................... 69 Appendix F: Simplification of double s ums without a multiplication factor ....................... 70 Appendix G: Simplification of double sums with a multiplication factor ............................ 71 Chapter 5. Appl ications of a Wertheim association activity coefficient model to methanol and ethanol - containing mixtures ................................ ................................ ................................ .......... 72 Introduction ................................ ................................ ................................ ........................ 72 Model theory and parameters ................................ ................................ ............................. 74 5.2.1. Association contribution ................................ ................................ ............................. 75 5.2.2. Combinatorial contribution ................................ ................................ ......................... 76 5.2.3. Residual contribution ................................ ................................ ................................ .. 78 5.2.4. EOS model parameters ................................ ................................ ............................... 82 5.2.5. Association parameters ................................ ................................ ............................... 82 Model comparisons ................................ ................................ ................................ ............ 85 5.3.1. Binary system with no association ................................ ................................ .............. 89 5.3.2. Binary systems with self - association ................................ ................................ .......... 92 5.3.3. Binary systems with self - and cross - associating systems ................................ ......... 103 5.3.4. Ternary systems ................................ ................................ ................................ ........ 110 5.3.5. General findings and recommendations ................................ ................................ .... 114 Association parameter value ................................ ................................ ............................ 117 Conclusions ................................ ................................ ................................ ...................... 120 APPENDIX ................................ ................................ ................................ ............................. 122 Chapter 6. Quantitative Analysis of Infr ared Spectra of n - Butanol + Cyclohexane with Quantum Chemical Calculations ................................ ................................ ................................ ................ 126 6.1. Introduction ................................ ................................ ................................ ...................... 126 6.1.1. Analysis of the O - H stretching band ................................ ................................ ......... 128 6.1.2. Motivation ................................ ................................ ................................ ................. 130 6.2. Computational methods ................................ ................................ ................................ ... 131 6.2.1. Molecular dynamics simulations ................................ ................................ .............. 131 6.2.2. Quantum mechanics simulations ................................ ................................ .............. 133 vii 6.3. Experimental methods ................................ ................................ ................................ ..... 134 6.4. Results and discussions ................................ ................................ ................................ .... 135 6.4.1. Processing and preliminary analysis of experimental IR spectra ............................. 135 6.4.2. Results from MD+QM ................................ ................................ .............................. 138 6.4.3. Scaling and further analysis of experimental IR spectra ................................ .......... 144 6.4.4. Universality of spectroscopic characteristics ................................ ............................ 148 6.5. Conclusions ................................ ................................ ................................ ...................... 150 Chapter 7. Integration of Quantum Calculatio ns and Spectroscopy for Wertheim Alcohol Association ................................ ................................ ................................ ................................ .. 151 7.1. Introduction ................................ ................................ ................................ ...................... 151 7.1.1. Hydrogen bonding and infrared spectroscopy ................................ .......................... 151 7.1.2. Thermodynamic modeling of hydrogen bonding ................................ ..................... 155 7.1.3. Motivation ................................ ................................ ................................ ................. 157 7.2. Experimental methods ................................ ................................ ................................ ..... 157 7.3. Results and discussions ................................ ................................ ................................ .... 158 7.3.1. IR band assignments ................................ ................................ ................................ . 164 7.3.2. IR peak fitting ................................ ................................ ................................ ........... 167 7.3.3. IR calculation of association parameters ................................ ................................ .. 175 7. 3.4. NMR calculations and results ................................ ................................ ................... 178 7.4. Summary of future work ................................ ................................ ................................ .. 183 APPENDIX ................................ ................................ ................................ ............................. 185 Chapter 8. Conclusions and Future Directions ................................ ................................ ........... 18 9 BIBLIOGRAPHY ................................ ................................ ................................ ....................... 193 viii LIST OF TABLES Table 3.1: Concentration based equilibrium constants ................................ ................................ . 30 Table 3.2: Equilibrium constant for primary alcohols ................................ ................................ .. 37 Table 3.3: three systems. ................................ ................................ ................................ ................................ 39 Table 4.1: Equations for sites in hypothetical system ................................ ............................ 52 Table 4.2: Key associations in the hypothetical system and their equilibrium constants ............. 53 Table 4.3: Other species formed in hypothetical system ................................ .............................. 54 Table 4.4: Association parameters and component densities ................................ ....................... 57 Table 4.5: Calculated monomer fraction values ................................ ................................ ........... 58 , = 0.2 and =0.6 ................................ ................................ ................................ ............. 59 Table 4.7 : Calculated parameters and intermediate values for activity coefficient calculations .. 65 Table 5.1: CPA parameters used to calculate ................................ ................................ ........ 83 Table 5.2: AspenPlus® databases used in calculations ................................ ................................ 85 Table 5.3: PC - SAFT pure component parameters ................................ ................................ ........ 85 Table 5.4: Data used for parameter regression ................................ ................................ ............. 86 Table 5.5: Total SSQ for systems studied in this work ................................ ................................ . 87 Table 5.6: Model performance in or der of best fit ................................ ................................ ........ 88 Table 5.7: Legend for plots in Model Capabilities section ................................ ........................... 89 Table 5.8: Association contribution to the total a ctivity coefficient at infinite dilution ............. 116 Table H.1: Regressed parameters values. ................................ ................................ ................... 123 Table H.2: Hayden - parameters. ................................ ................................ ........... 125 ix Table 6.1: MD simulation details ................................ ................................ ................................ 132 Table 6.2: Number of each species and bond analyzed with QM calculat ions .......................... 138 Table 6.3: Average calculated vibrational frequencies and integrated absorption coefficients . 140 Table 6.4: Areas and sta ndard deviations for the absorbance and scaled absorbance bands ...... 147 Table 7.1: Frequency range defined as hydroxyl band for each spectrum ................................ . 158 Table 7.2: Number of each species and bond analyzed with QM calculations for the distribution plot in Figure 7.7 and Figure 7.8 ................................ ................................ ................................ . 164 Table 7.3: Concentration of free ends ( ) in mol/mL calculated from fitted Gaussian peaks ................................ ................................ ................................ ................................ ..................... 175 Table 7.4: Regressed parameters from NMR analysis and corresponding CPA parameters for comparison ................................ ................................ ................................ ................................ .. 182 Table 8.1: Accomplished and future goals for research project ................................ ................. 190 x LIST OF FIGURES Figure 1.1: NRTL representat ion of VLE and LLE of methanol (1) cyclohexane (2). ............... 2 Figure 1.2: Vapor pressure as a function of temperature for ethanol calculated using two sets of CPA parameters. ................................ ................................ ................................ ............................. 5 Figure 3.1: 2B association scheme. ................................ ................................ ............................... 27 Figure 4.1: 2B association scheme. ................................ ................................ ............................... 49 Figure 4.2: H ypothetical molecules for analysis. ................................ ................................ .......... 50 Figure 4.3: Bonding schemes for hydrogen bonding in a mixture of E , F and G . ........................ 53 Figure 4.4: Counting species in chemical theory system ................................ .............................. 55 Figure 4.5: True monomer fraction of components in system of study and fractions of acceptors and donors ................................ ................................ ................................ ................................ ..... 60 Figure 4.6: Distribution of species with component E when =0.2 ................................ ........... 60 Figure 4.7: Distribution of species with component F when =0.2 ................................ ........... 61 Figure 4.8: Distribution of species with component G when =0.2 ................................ ........... 61 Figure 4.9: True mole fraction of methanol chains calculated with 4 values of ............... 64 Figure 4.10: Calculated values of the activity coefficient contributions ................................ ...... 66 Figure 5.1: Phase equilibria diag ram for n - heptane + cyclohexane at 101.3 kPa ......................... 90 Figure 5.2: Heat of mixing for n - heptane + cyclohexane at T=298.14 K and 101.3 kPa ............. 90 Figure 5.3: Calculated and experimental activity coefficients for n - heptane + cyclohexane at 101.3 kPa ................................ ................................ ................................ ................................ ................. 91 Figure 5.4: Percent errors of K - ratios for n - heptane + cyclohexane system ................................ 91 Figure 5.5: Phase equilibria diagram for methanol + n - heptane at 101.3 kPa .............................. 92 Figure 5.6: Heat of mixing for methanol + n - heptane at T=2 98.15 K and 101.3 kPa .................. 93 Figure 5.7: Calculated and experimental [4] activity coefficients for methanol + n - heptane at 101.3 kPa ................................ ................................ ................................ ................................ ................. 93 xi Figure 5.8: Percent errors of K - ratios for methanol + n - heptane system ................................ ...... 94 Figure 5.9: Phase equilibria diagram for methanol + cyclohexane at 101.3 kPa .......................... 95 Figure 5.10: Heat of mixing for methanol + cyclohexane at T=416.29 K and 1900 kPa ............. 95 Figure 5.11: Calculated and experimental activity coefficients for methanol + cycl ohexane at 101.3 kPa ................................ ................................ ................................ ................................ ................. 96 Figure 5.12: Percent errors of K - ratios for methanol + cyclohexane system ............................... 96 Figure 5.13: Phase equilibria diagram for methanol + n - pentane at T=397.7 K .......................... 97 Figure 5.14: Calculated and experimental [145] activity coefficients for methanol + n - pentane at T=397.7 K ................................ ................................ ................................ ................................ ..... 97 Figure 5.15: Percent errors of K - ratios for methanol + n - pentane system ................................ .... 98 Figure 5.16: Phase equilibria diagram for ethanol + cyclohexane at 101.3 kPa ........................... 99 Figure 5.17: Heat of mixing for ethanol + cyclohexane at T=298.14 K and 101.3 kPa ............... 99 Figure 5.18: Calculated and experimental activity coe fficients for ethanol + cyclohexane at 101.3 kPa ................................ ................................ ................................ ................................ ............... 100 Figure 5.19: Percent errors of K - ratios for ethanol + cyclohexane system ................................ 100 Fi gure 5.20: Phase equilibria diagram for ethanol + n - decane at T=338.17 K ........................... 101 Figure 5.21: Phase equilibria diagram for ethanol + n - decane at 101.3 kPa .............................. 101 Figure 5.22: Heat of mixing for ethanol + n - decane at T=298.15 K and 101.3 kPa .................. 102 Figure 5.23: Phase equilibria diagram for methanol + water at 101.3 kPa ................................ . 104 Figure 5.24: Phase equilibria diagram for methanol + water at T=298.14 K (bottom) and T=308.14 K (top) ................................ ................................ ................................ ................................ ......... 104 Figure 5.25: Calculated and expe rimental [146,147] activity coefficients for methanol + water at T=298.14 K, T=308.14 K and P=101.3 kPa from top to bottom ................................ ................ 105 Figure 5.26: Percent errors of K - ratios for methanol + water syste m ................................ ........ 106 Figure 5.27: Phase equilibria diagram for methanol + ethanol at T=298.14 K (bottom) and T=413.13 K (top) ................................ ................................ ................................ ........................ 107 Figure 5.28: Ph ase equilibria diagram for methanol + ethanol at T=298.14 K .......................... 108 xii Figure 5.29: Calculated and experimental activity coefficients for methanol + ethanol at T=298.14 K (bottom) and T=413.13 K (top) ................................ ................................ .............................. 108 Figure 5.30: Percent errors of K - ratios for methanol + ethanol system ................................ ...... 109 Figure 5.31: Phase equilibria diagram for methanol + n - heptane + cyclohexane at 101.3 kPa and 298.15 K ................................ ................................ ................................ ................................ ...... 110 Figure 5.32: K - ratios for methanol (a) + n - heptane (b) + cyclohexane (c) system at 101.3 kPa and 298.15 K ................................ ................................ ................................ ................................ ...... 111 Figure 5.33: Phase equilibria diagram for methanol + ethanol + cyclohexane at 101.3 kPa and 298 K ................................ ................................ ................................ ................................ .................. 112 Figure 5.34: K - ratios for methanol (a) + ethanol (b) + cyclohexane (c) system at 101.3 kPa and 298 K ................................ ................................ ................................ ................................ ........... 113 Figure 5.35: Comparison of phase equilibria fits with two sets of CPA parameters for ethanol + cyclohexane at 101.3 kPa ................................ ................................ ................................ ............ 118 Figure 5.36: Comparison of phase equilibria fits with two sets of CPA parameters for ethanol + n - decane at T=338.17 K ................................ ................................ ................................ ................. 118 Figure 5.37: for the self - association of the primary alcohols at 298.15 K calculated by fitting to different properties. ................................ ................................ ................................ ................. 119 Figure 6.1: Types of covalent O - H bonds ................................ ................................ ................... 130 Figure 6.2: Bond distribution in trimers calculated with MD simulations for an equimolar n - butanol + cyclohexane mixture ................................ ................................ ................................ ... 133 Figure 6.3: Experimental O - H IR ba nd for a 0.1 mole fraction of n - butanol in cyclohexane .... 136 Figure 6.4: Total absorbance band area of the O - H band for n - butanol + cyclohexane data as a function of the mole fraction for 4 compo sitions ................................ ................................ ........ 137 Figure 6.5: Hydroxyl stretching frequencies and integrated absorption coefficients for two concentrations of n - butanol + cyclohexane mixtures calculated from QM simulations ............ 141 Figure 6.6: Hydroxyl stretching frequencies and integrated absorption coefficients for other systems calculated from QM simulations ................................ ................................ ................... 144 Figure 6.7: Scaled absorbance spectra for n - butanol + cyclohexane at = 0.1 ................. 145 Figure 6.8: Scaled absorbance band area of the O - H band for n - butanol + cyclohexane data as a func tion of the mole fraction of 4 experimental compositions ................................ ................... 146 xiii Figure 6.9: Relationship between the covalent O - H bond length and vibrational frequency calculated for equimolar ethanol + cyclohexa ne mixture ................................ ........................... 148 Figure 6.10: Comparison of calculated NMR and IR characteristics for equimolar ethanol + cyclohexane mixture ................................ ................................ ................................ ................... 149 Fi gure 7.1: Types of covalent OH bonds ................................ ................................ .................... 153 Figure 7.2: Hydroxyl stretching frequencies and integrated absorption coefficients for n - butanol + and ethanol + cyclohexane mixtures calculated from QM simulations ................................ ...... 154 Figure 7.3: Raw (top) and scaled (bottom) O - H IR band for a 0.0469 mole fraction of n - butanol in cyclohexane ................................ ................................ ................................ ................................ . 160 Figure 7.4: Raw (top) and scaled (bottom) O - H IR band for a 0.0678 mole fraction of n - butanol in cyclohexane ................................ ................................ ................................ ................................ . 161 Figure 7.5: Raw (top) and scaled (bottom) O - H IR band for a 0.0817 mole fraction of n - butanol in cyclohexane ................................ ................................ ................................ ................................ . 162 Figure 7.6: Raw (top) and scaled (bottom) O - H IR band for a 0.100 mole fraction of n - butanol in cyclohexane ................................ ................................ ................................ ................................ . 163 Figure 7.7: Smoothed normalized distributions of bond types from QM simulations of n - butanol + and ethanol + cyclohexane mixtures. ................................ ................................ .......................... 165 Figure 7.8: Smoothed normalized distributio ns of bond types from QM simulations with from dimers separated from others. ................................ ................................ ................................ ..... 166 Figure 7.9: Peaks fitted to scaled IR spectra for n - butanol + cyclohexane mixture at =0.0469. ................................ ................................ ................................ ................................ ..................... 168 Figure 7.10: Peaks fitted to scaled IR spectra for n - butanol + cyclohexane mixture at =0.0678. ................................ ................................ ................................ ............................ 169 Figure 7.11: Peaks fitted to scaled IR spectra for n - butan ol + cyclohexane mixture at =0.0817. ................................ ................................ ................................ ............................ 170 Figure 7.12: Peaks fitted to scaled IR spectra for n - butanol + cyclohexane mixture at =0.100. ................................ ................................ ................................ ................................ ..................... 1 71 Figure 7.13: Scaled band area as a function of apparent n - butanol concentration ..................... 172 Figure 7.14: Bond distributions calculated from IR peak fitting ................................ ................ 173 Figure 7.15: Concentration of bonds vs. mole fraction of n - butanol ................................ ...... 174 Figure 7.16: Calculated and fitted concentrations of free end ( hydroxyls ..................... 176 xiv Fi gure 7.17: Derived association parameter values from IR analysis (solid line) compared to CPA values (dashed line) for the self - association of n - butanol ................................ .................. 177 Figure 7.18: O - H proton chemic al shift as a function of n - butanol mole fraction for n - butanol + cyclohexane mixtures at 26.3 C. ................................ ................................ ............................... 180 Figure 7.19: O - H proton chemical shift as a function of n - butanol mole fraction for n - buta nol + cyclohexane mixtures at 34 C. ................................ ................................ ................................ .. 181 Figure 7.20: O - H proton chemical shift as a function of n - butanol mole fraction for n - butanol + cyclohexane mixtures at 41 C. ................................ ................................ ................................ .. 181 xv KEY TO SYMBOLS AND ABBREVIATIONS a interaction parameter for molecules in a mixture, calculated with mixing rules interaction parameter for i number of molecules of a single component in Chapter 3; activity of species i otherwise interaction parameter for interactions between molecules of components i and j CPA parameter for attractive term fitted spectroscopic parameter fitted residual term parameter A molar Helmholtz energy (extensive if accompani ed by underbar); absorbance in Chapters 6, 7 absorbance of species i scaled absorbance fitted residual term parameter AVEC association with variable equilibrium constant size parameter of mixture, calculated with mixing rules size parameter for i number of molecules of a single component fitted spectroscopic parameter fitted residual term parameter B free free contribution to the second virial coefficient Intermediate CPA parameter for attractive t erm fitted spectroscopic parameter concentration of component i (apparent unless accompanied by subscript T) apparent dimensionless concentration of component i xvi dimensionless true monomer concentrations of componen t i in a mixture dimensionless true monomer concentrations of component i in a pure solution of component i CLAM continuous linear association model CPA cubic - plus - association model fitted spectroscopic parameter DFT density functi onal theory EOS equation of state ECR ESD Elliott - Suresh - Donohue model constant in Eq. ( 3 . 14 ) fugacity of component i fitted spectroscopic parameter function defined by Eq. ( 3 . 18 ) in Chapter 3; radial distribution function otherwise fitted spectroscopic parameter molar Gibbs energy (extensive if accompanied by underbar) function defined by Eq. ( C. 1 ) fitted spectroscopic parameter molar enthalpy (extensive if accompanied by underbar) IR infrared spectroscopy reaction rate constant for forward reaction reaction rate constant for reverse reaction xvii binary inte raction parameter K equilibrium constant fitted association parameter with units of volume KLL K - ratio in the liquid - liquid region KVL K - ratio in the vapor - liquid region pathlength generic notation for a regressed property length of the hydroxyl bond LLE liquid - liquid equilibrium LACT linear association with cyclic trimer LCST lower critical solution temperature number of segments per chain in PC - SAFT MD molecular dynamics molecular weight of species i num ber of moles number of non - bonded sites of type total number of sites of type in solution n* effective number of carbon atoms in an ether homomorph molecule number of molecules of compo nent i local composition, number of component i molecules around molecule j number of sites identical to site i of type B on the site host synonymous to , introduced for disambiguation in Eq. ( 4 . 1 ) NMR nuclear magnetic resonance spectroscopy xviii NPT condition holding number of moles, pressure and temperature constant NRTL non - random two liquid model NVT condition holding number of moles, volume and temperature constant P probability in Chapter 4, pressure otherwise Gaussian peak fitted to IR vapor pressure of component i PC - SAFT perturbed chain statistical associating fluid theory PM1, PM2 CPA pure component parameters from Kontogeorgis et al . [1] and [2] respectively PME Particle Mesh Ewald QM quantum mechanics r i ratio of reference molar density, ( taken to be that of methanol at 303.15 K) to the molar density of component i, R universal gas constant standard deviation molar entropy of solution molar entropy of component i SAFT statistical associating flui d theory model SG Staverman - Guggenheim model SH Scatchard - Hildebrand model SSQ sum of squares error t time T temperature molar internal energy (extensive if accompanied by underbar) xix UNIFAC UNIQUAC functional - group activity coefficient model UNIQUAC univ ersal quasichemical model V i molar volume of pure component i V molar volume of solution, extensive volume if accompanied by underbar VLE vapor - liquid equilibrium VLLE vapor - liquid - liquid equilibrium WAG Wertheim association gamma model apparent l iquid mole fraction of component i mole fraction of component i molecules around a molecule of component j fraction of sites identical to site i of type B that remain non - bonded at equilibrium vapor phase mole fraction of s pecies i true mole fraction (for monomer in Chapter 4) Z compressibility factor Greek Symbols type of O - H bond in Chapters 6 and 7 given in Figure 6 . 1 non - randomness parameter in NRTL model function defined by Eq. ( E. 2 ) type of O - H bond in Chapters 6 and 7 given in Figure 6 . 1 association volume parameter activity coefficient; type of O - H bond in Chapters 6 and 7 given in Figure 6 . 1 type of O - H bond in Chapters 6 and 7 given in Figure 6 . 1 solubility parameter in Scatchard - Hildebrand Werth xx displacement of hydrogen atom in x - direction displacement of hydrogen atom in y - direction displacement of hydrogen atom in z - direction observed NMR chemical shift relative to monomer NMR chemic al shift NMR chemical shift of hydrogen bonded proton relative to monomer NMR chemical shift depth of pair potential in PC - SAFT fitted association interaction energy parameter molar absorption coefficient molar absorption coefficient of species i type of O - H bond in Chapters 6 and 7 given in Figure 6 . 1 packing fraction; type of O - H bond in Chapters 6 and 7 given in Figure 6 . 1 thermal correction factor effective association volume in PC - SAFT rearrangement of ( 5 . 11 ) stoichiometric number of species i w avenumber, = frequency ( )/speed of light ( NMR chemical shift molar solution density (apparent unless accompanied by subscript T) molar density of component i (apparent unless accompanied by subscript T) interaction energy parameter between components i and j size parameter, i.e. contact distance temperature - independent segment diameter in PC - SAFT xxi fugacity coefficient of pure component i volume fraction of component or species i (apparent unless accompanied by subscript T) modified apparent volume fraction of component or species i integrated absorption coefficient i true dimensionless concentration of al l chains that end in i local composition parameter defined by Eq. ( 5 . 7 ) Subscripts 0 pertaining to apparent species activity - coefficient a activity - based alc pertaining to the alcohol avg averaged C concentration - based CPA pertaining to CPA Crit critical property D dimer f pertaining to the formation of a reaction product h host pertaining to component hosting the site at maximum peak height p ertaining to site j that can bond to mix mixture property xxii M monomer constant number of moles for all components except component Nmer complex (hydrogen - bonded species) p constant power P at constant pressure T pertaining to true spec ies or at constant temperature Tri trimer at constant volume x mole fraction based Superscripts standard state property 0 property for pure component assoc att a ttractive contribution interaction between site i of type B and site j of type E c cyclic species chem chemical theory comb combinatorial contribution exp experimental value E excess property Flory heory hard sphere xxiii ig ideal gas is ideal solution l linear species L pertaining to a liquid Mod Flory Nagata NRTL calculated with NRTL constant for Eq. ( 5 . 6 ) ref pertaining to a reference state rep repulsive contribution res residual contribution property at saturation SH calculated with Scatchard - Hildebrand T total pertaining to a vapor Wilson Special Notat ions as in extensive property as in partial molar property as in property of component in a mixture [ i ] apparent concentration of component i 1 Chapter 1. Introduction to T hermodynamic M odeling of A ssociating S ystems Thermodynamic models play an essential role in the chemical industry. These tools provide valuable information on the qualitative and quantitative responses of a system to variations in temperature, pressure and composition . As such, they are necessary in the design of process equi cannot be investigated with bench - scale experiments. Of particular importance is the use of thermodynamic models for the design and operation of separation units, such as distillation the optimization of each component of a separation process , including the accuracy and robustness of the thermodynamic models, is an active are a of research. Dow Chemical and Fluor manufacture large quantities of alcohols, ketones and esters. Honeywell UOP is interested in biorefining pyrolysis oils to fuels, a process that includes three - phase mixtures that are difficult to model. In cases like these where one or more components in a system can hydrogen bond, the relatively strong associations within and between molecules complicate thermodynamic model ing. Indeed, traditional models do not specifically account for hydrogen bonding between molecules, but instead lump effects into dispersion terms of the equations. For example, the non - random two liquid (NRTL) [3] model is commonly used in industry due to its relative simplicity and flexibility. This model represents molecular asymmetry with up t o three parameters per binary, two for pair interactions at infinite dilution and one to adjust model to deviations from ideality at intermediate compositions. However, because the binary parameters are intended to represent relatively nondirectional dispe rsion interactions the model can fail to accurately represent components that hydrogen bond. This means that, when extrapolated to new 2 sets of conditions, NRTL can perform poorly and can therefore not be used with confidence . This is illustrated in Figure 1 . 1 for the vapor - liquid (VLE) and liquid - liquid equilibria (LLE) of methanol in cyclohexane. When the VLE envelope is fitted (solid lines), the same parameters yield a poor prediction of the LLE, incorrectly indic ating a three - phase vapor - liquid - liquid (VLLE) condition which would occur where the LLE binodal lines extrapolate to the VLE. Conversely, when the LLE is fitted (dashed lines), the model falsely predicts a higher temperature azeotrope and improperly shift s the VLE. Figure 1 . 1 : NRTL representation of VLE and LLE of methanol (1) cyclohexane (2) . Solid lines are fitted to the VLE envelope ( , , and ) and dashe d lines are fitted to the LLE envelope ( , , and ) . Markers are experimental VLE [4] and LLE [5] phase equilibria dat a. To make models more flexible and robust, it has been demonstrated [6,7] that hydrogen bonding is best calculated explicitly as a separate term. 270 280 290 300 310 320 330 340 350 360 0 0.2 0.4 0.6 0.8 1 Temperature (K) Mole fraction of methanol VLE VLE LLE 3 1.1. Background Hydrogen bonding is a type of dipole - dipole interaction that occurs between a hydrogen atom bonded to a 'p roton donor' and a neighboring strongly electronegative 'proton acceptor'. While weak compared to a covalent bond , this interaction is responsible for the high er melting and boiling points of water and alcohols relative to hydrocarbons with similar molecul ar weights. - - to hydrogen bonding between two molecules of the sa me component, such as between two methanol molecules or between two water molecules. Cross - association describes hydrogen bonding between molecules of two different components which are also able to self - associate, such as an interaction between a water mo lecule and a methanol molecule. The term solvation is used in cases where one or both components involved in a hydrogen bond do not self - associate. Examples of solvating systems are water + ethers or chloroform + ethers. Though these definitions describe d istinct behavior, the terms are often used loosely in literature. Frequently, the term collectively and we adopt this convention in the present document. To provid e context for current state of the art in association modeling, we describe two categories under which most thermodynamic models fall: equations of state (EOS) and activity coefficient ( models. Thermodynamic modeling is based on quantifying departures or difference from idealized mixture states. The primary difference between the two classes is the choice of idealized state. Ideality, as it refers to ideal gases or solutions, is a simplified hypothetical state which is used as a conceptual reference poi nt in thermodynamic model design. From this starting point, 4 contributions are incorporated to calculate the differences between ideal and real solution dynamics. In many cases, the major contribution to non - ideality is due to hydrogen bonding. 1.1.1. Equations o f state (EOS) Equations of state (EOSs) use departures from ideal gases to model mixture behavior. An ideal gas is a collection of point particles which do not interact with one another by collision, repulsion or attraction. The fugacity of each component of a gas mixture, , is given by [8] 1 . 1 where is the vapor phase mole fraction of component i , is the total pressure. The component fugacity coefficient, , quantifies the departure from an ideal gas mixture, and is unity for an ideal gas. Equations of state (EOS), such as the van der Waals and P eng - Robinson equations of state, calculate deviations from this idealized state by adding the effects of repulsive and attractive forces and association in solution. As such, the fugacity coefficient is separated as 1 . 2 where is the compressibility factor and can also be separated into the same contributions. 1 . 3 Classical EOSs omit the association term entirely. However, newer models such as the S tatistica l A ssociating Fl uid T heory (SAFT) [9] , Elliott - Suresh - Donohue (ESD) [10] and C ubic - P lus - A ssociation (C PA) [11] were developed to include hydrogen bonding explicitly. While the theoretical basis of these models is solid, and they are popular in academics, there are th ree barriers that impede their wide adoption in the chemical industry. The first barrier that limits wide spread use is that the parameters for these models are determined by fitting to a large quantity of pure component vapor pressure and density experim ental data, 5 which are frequently unavailable to industry over a suitable span. When new compounds are developed, data is expensive to obtain over wide ranges, and limited data results in significant uncertainty in parameter fitting. The second barrier to wider adoption is that dispersion forces and the association forces are both demonstrated for ethanol in Figure 1 . 2 using calculations of the pure component vapor pressure as a function of temperature. Figure 1 . 2 : Vapor pressure as a function of temperature for ethanol calculated using two sets of CPA parameters . Blue markers are experimental data, black dotted line are calculated with parameters from Kontogeorgis (2006) [1] (referred to in Chapter 5 as PM1) and dashed lines are calculated with parameters from Kontogeorgis (2010) [2] (referred to in Chapter 5 as PM2). The bottom plot is a zoom - in of the top plot. -2 0 2 4 6 8 0.0018 0.0024 0.003 0.0036 0.0042 0.0048 0.0054 log 10 (P sat (Pa)) 1/T (K - 1 ) 3.5 4 4.5 5 5.5 6 0.0024 0.0026 0.0028 0.003 0.0032 0.0034 log 10 (P sat (Pa)) 1/T (K - 1 ) Exp Data Prms from Kontogeorgis 2006 Prms from Kontogeorgis 2010 6 The two black lines are both calculated using the CPA EOS and the only differe nce is the fitting method used to obtain the pure component parameters. For the dotted line, the parameters were fitted to pure vapor pressure and density (PM1) [1] and, for the dashed line, parameters were fitted to vapor pressure, density, and spectroscopic data (PM2) [2] . The numerical values of the parameters in PM2 are up to 37% different than those in PM1. In both fits presente d in Figure 1 . 2 , the CPA parameters were fitted to data over a wide temperature range. Despite this, the PM1 parameters clearly fit the vapor pressure much better than PM2 parameters. Indeed, the PM2 fit underesti mates the vapor pressure of ethanol by 22 kPa at T=338.17 K when the constraint is added to also represent spectroscopic data. The challenge in fitting indicates that several reasonable solutions exist for the fitting problem and that it is under - specified , and still may not represent all properties accurately. Even when fitted to a large amount of data, the parameters are sometimes not known with confidence. Making small adjustments to one parameter requires refitting of all parameters. The third obstacle for more wide - spread adoption of EOSs is that the ability to represent mixtures is sensitive to accuracy of fitting the pure component property data. For example, if the parameters fit vapor pressure poorly at purity such as ethanol with the PM2 parameter set, the resulting mixture phase diagram will be skewed and fail to match experiment, even if the mixture deviations from mixture ideality are correct. This difficulty is further discussed in Chapter 5. 1.1.2. Activity coefficient models In contrast to EOSs, ac tivity coefficient ( ) models, such as NRTL and UNIversal QUAsiChemical (UNIQUAC) [12] , calculate deviations from an ideal solution state and are commonly used for liquid states. For this approach, component fugacities of each component i in solution are calculated relative a standard state liquid fugacity, , often set to be that of a pure component. 7 1 . 4 where is the mole fraction of component i . The ideal solution model requires that the activity coefficien t, , is unity. More specifically, the ideal solution model requires several important constraints, including: 1. No excess volume , where and are the mixture molar volu me and the component molar volume respectively , an d no excess entropy, where and are the mixture molar entropy and the component molar entropy respectively . These constraints are approached when m olecules in a mixture have similar size s and sh ape s . 2. No excess energy or This constraint requires that t he energy of all intermolecular interactions, including hydrogen bonding, are the same regardless of the molecule types involved across the composition range. This constraint is approached when mixing methanol and ethanol for example. Deviations from these ideality constraints are corrected in activity coefficient models by combinatorial, residual and association contribution respectively. Models are typically developed for the logarithm s of the contributions, which are added according to [8] : 1 . 5 Compared to EOSs, activity coefficient models are not limited by the first and third challenges outlined in S ection 1.1.1 . First, models can be used with limited pure component data because vapor pressure and density data are required only for the limited range where the models must be applied. Secondly, pure component fugacities and vapor pressure s can be represented by correlations separate from the mixture model, decoupling the fitting of the pure component and mixture properties. 8 Despite these advantages, activity coefficient models that capture effects of hydrogen bonding have not gained wide acceptance. Until the mid - 1990s , models primarily used a chemical theory approach where the association is represented by chemical species [13 16] . Chemical theory is difficult to generalize, which creates a barrier for wide - spread use. Fu and Sandler [6] develop expressions that are readily applied in a generalized way. Development of activity mo dels has languished while most development effort has focused on equations of state. 1.1.3. Association modeling Association perturbation theory, which is the method used in CPA and th e SAFT families of equations. While such similarity is that both methods have two key quantities: an association constant and another variable that is related - bonded at equilibrium. A common challenge between EOS and activity coefficient models is discerning the value for these variables. In fact, published association constants can vary up to a factor of four to five depending on the source. A fundamentally satisfying method to obtain accurate association constants is to quantify the real monomer fraction in solution and fit the association parameters to the monomer fraction. Some efforts have been made to determine monomer fractions from molecular dynamics simulations using atomistic force fields. While more rigorous than fitting to macroscopic density and vapor pressure data, the monomer fractions and hydrogen bonding depend significantly on the force field employed. Kwac and Geva [17] show that the relative population of monomers can vary by up to 13% between various force fields. A better approach is to probe the microscopic 9 nature of the interactions with spectroscopic tool s such as infrared (IR) or nuclear magnetic resonance (NMR) spectroscopy. The advantage of these techniques lies in the measurement of the experimental response due to the types and relative quantities of microscopic species existing in solution. However, attainment of high - quality information from these techniques remains challenging. For example, this work will show that the IR extinction coefficient varies by approximately 20 times for different types of hydroxyl groups in solution. Due to overlapping pe aks, individual species have not previously been determined experimentally with confidence. NMR provides a time - averaged measure of hydrogen bonding and, in cases with multiple equilibrating species, the measurements provide only broad insight into the spe cific distribution of hydrogen bonded species in solution. 1.2. Research o bjectives The project developed in this dissertation integrates multiple tools to improve the thermodynamic modeling of associating systems: a new activity coefficient model, spectroscopy and MD + QM simulations. The complimentary relationship between these facets provides new quantitative insights that overcome the limitations and exceed the capabilities of the individual tools acting alone. First, a flexible thermodynamic activity coeffi cient model is developed that incorporates from IR and NMR spectra. Finally, to overcome the challenges associated with spectral analysis, MD + QM simulations are used to estimate the extinction coefficients of various species and predict their distribution in solution. Ultimately, the objective of this work is to advance knowledge of hydrogen bonding by quantifying its effects at a molecular level to lend physical significance to engineering models for industry. 10 1.3. Dissertation o utline This remainder of this document is organized as follows: Chapter 2 describes the relation between the concentration - based equilibrium constant, and the activity - based equilibri um constant, . Equilibrium reactions and chemical theory are often described by which is assumed to be independent of composition in the liquid phase. However, from the thermodynamic perspective, it is actually that remains constant with composition as it is directly related to the Gibbs energy. The chapter explores the relationship between the two types of equilibrium constants, and their use in reaction engineering. In Chapters 3 and 4, the two most prevalent association theories in numerically and analytically identical under reasonable approximations. We demonstrate a direct relationship between the concentration - based equilibrium constant used in chemical theory, and the Wertheim association constant by relating them both to the activity - based equilibrium constant, . Moreover, an activity coefficient association contribution ( ) model, named th e Wertheim Association Gamma (WAG), is developed to provide generalized method to integrate association into thermodynamic models. Chapter 3 considers associating species in an inert solvent. In Chapter 4, the analysis is extended to systems that solva te in addition to forming chains. Existing literature in this area has rarely detailed the species in chemical theory. Therefore, objectives of Chapters 3 and 4 are to elucidate how various species must be counted in chemical theory and to demonstrate the relationship with Wertheim theory. The WAG model developed in Chapter 3 is used in Chapter 5 to model the behavior of methanol and ethanol - containing binary and ternary systems of industrial significance. The purpose of this work is to demonstrate the capa bilities and limitations of the functional form of the when 11 combined with traditional physical models such as NRTL [3] , Wilson [18] , and Scatchard - Hildebrand [19,20] and compare with the SAFT and CPA models. To this end, the values o f the association constants ( are set to values determined from published CPA parameters. In Chapters 6 and 7, fundamental studies are leveraged to obtain more physically significant values of . First, a successive mol ecular dynamics and quantum mechanics (MD + QM) approach to investigating molecular association is outlined. Molecules are classified according to their participation in hydrogen bonding, and their response to excitation by infrared is studied. A universal map is created that relates the hydroxyl IR band characteristics and the bond length of the hydroxyl in alcohol monomers and clusters. Then, the QM maps are used to analyze experimental IR spectra to calculate for butanol + cyclohexane. In the final Chapter 8, the main conclusions of the project are summarized and recommendations for future work are outlined. 12 Chapter 2. Relation of Concentration - based Equilibrium Constants to Activity - based Equilibrium Constants in the Liquid Phase Concentration - based equilibrium constants are used widely in modeling chemical kinetics for elementary reactions in both the gas phase and liquid phase. However, activity - based equilibrium constants are the rigorous values ba sed on thermodynamic minimization of the Gibbs energy. The relation between the concentration - based constants and the activity - based constants is well - documented in textbooks for the ideal gas phase. This research note shows the relation between the equili brium constants in the liquid phase and discusses conditions under which the concentration - based equilibrium constant is independent of composition. 2.1. Introduction Characterization of chemical equilibria is necessary for chemical reactor design, and of inter est for understanding limitations of chemical conversion. A requirement for chemical equilibrium is that the Gibbs free energy of the reaction is minimized, and the standard state Gibbs energy change is related to the activity - based equilibrium constant, [8] . As such, from a fundamental thermodynamic perspective, is independent of composition for a given reaction. However , in practice, this compositional independence is commonly assumed for the concentration - based equilibrium constant, , for reactions in both the gaseous and liquid phases. T he use of concentrations in the liquid phase is often explained using the em pirical law of mass action rather than on a basis of thermodynamic principles [21 23] . In this work, we show why this approach is valid for ideal gases and explore the conditions under which the approximation is reasonable for liquids. 13 2.2. Gibbs energy and the activity - based equilibrium constant For the purposes of this communication, we consider the reaction components E and F combining into a product composed of and . 2 . 1 Commonly, reactions are generalized using stoichiometric coefficients, and stoichiometric numbers, [21] . Thus, the stoichiometric coefficient s for this reaction are e, f, and 1 respectively for the three species, and the stoichiometric numbers are , , and . While there is only one product for this simple example, the results set forth here are not limited to single - product reactions. Thermodynamic analysis focuses on the activity - based equilibrium constan t defined by 2 . 2 where is the activity of species i and depends on the standard state. The activity - based equilibrium constant is dimensionless, independent of composition and depends only on temperature. The equilibrium constant value is calculated via 2 . 3 where is the standard Gibbs energy of r eaction and is the standard state Gibbs energy of formation of component i . The value of the equilibrium constant depends on the standard state used for each component and the value will be different when ideal gas standards states are used compared to when liquid standards states are used. Commonly, all the components use the same standard state, though this is not required. In many cases, pure thermodynamic quantities are not 14 known, and the value of the equilibrium constant is determined em pirically from applying Eq. ( 2 . 2 ) or approximations. The activity in Eq. ( 2 . 2 ) is given by 2 . 4 where is the component fugacity and is the standard state fugacity. A standard state is specified at a particular composition, state of aggregation, and typically chosen to be at the temperature of the reacting system and a specified standard state pressure, . 2.3. Gas phase a ctivity - b ased e quilibrium c onstant When the reactants and products are in the gas phase, the fugacity is calculated using an equatio n of state, , where is vapor mole fraction and is the component fugacity coefficient. Using the common standard state of a pure ideal gas at a pressure of 1 bar, bar, and for the reaction of Eq. ( 2 . 1 ), inserting activities from Eq. ( 2 . 4 ) into Eq. ( 2 . 2 ) gives 2 . 5 w here commonly is omitted from the expression and is expressed in bar to assure that remains dimensionless. When the gas phase is ideal, t he fugacity simplifies by neglecting gas phase nonidealities, , and recognizing the molar concentration for species i , and denoting it by square brackets , : 2 . 6 The expression in Eq. ( 2 . 6 ) is only valid for ideal gases as noted with the superscript ig . The concentration - based equilibrium constant found within the equation is defined as 2 . 7 15 This relation is s hown widely in many texts and developed only for ideal gases and/ or used without proof for condensed phases [8,21 25] . Si nce is independent of composition, it is obvious from Eq. ( 2 . 6 ) that is independent of concentration for ideal gas species. However, for liquid phases, the ideal gas law does not apply and thus, . Therefore, is not the correct conversion factor between the liquid phase and values and the general relation for condensed phases is not obvious. The following section demonstrates the relation between and in condensed phases and explains the conditions under which liquid - phase can be expected to be independent of composition. 2.4. Liquid p hase a ctivity - b ased e quilibrium c onstant One method to explore the composition dependence of in the liquid phase is to use a standard state of 1 M in Eq. ( 2 . 3 ), and to use the concentration scale for activities in Eq. ( 2 . 4 ). However, activity relations for nonele ctrolytes on the concentration scale are uncommon. This work will use the mole fraction scale to understand the relation between and and the conditions where may be expected to be independent of composition in the liquid phase. The use of also arises from reaction kinetics. Chemical kinetics are routinely modeled using concentrations to represent the d riving force [21 23] for the reaction rate. For the reaction in Eq. ( 2 . 1 ), the commonly - used concentration - based rate expression for an elementary reversible reaction in a batch reactor is 2 . 8 where and are the forward and reverse rate constants, and are the rates of change of E and concentrations, respectively. The term is the rate of change for concentration of . Rate laws represent the non - equilibrium reaction phenomena, and a lengthy 16 review regarding the use of the law of mass action is av ailable [26] . For the purposes of this note, we accept the form of Eq. ( 2 . 8 ) as suffic ient for engineering application for ideal gases. If valid, since the forward and reverse rates are equal at equilibrium, then 2 . 9 In the liquid phase, using a liquid standard state, the activity is where and are evaluated according to and at P and P o respectively. Ignoring the minor effect of pressure on liquid fug acity, the equilibrium constant for the example reaction of Eq. ( 2 . 1 ) results in [8] 2 . 10 To introduce concentrations for condensed phases, molar densities are typically used according to 2 . 11 where is the mixture molar density. Thus, the relation between and is 2 . 12 where 2 . 13 Because is independent of composition, and molar density and activity coefficients change with composi tion, then may be expected to display composition dependence via Eq. ( 2 . 12 ), which is in conflict with the frequent empirical use of for modeling of reactions that approach reaction equilibrium . Cer tainly, activity - based kinetics can model the composition dependence [27,28] , but improved clarity of the relation between the and would benefit the kinetics and thermodynamics communities. 17 Within chemical engineering, activity coefficients are commonly modeled by [12,29] 2 . 14 where comb [30] or some variant, and res indicates the energetic and other residual ef fects that are commonly modeled with adjustable parameters [12,29] . Note that Eq. ( 2 . 14 ) combined with Eq. ( 2 . 10 ) becomes 2 . 15 where the superscript L indicates liquid phase and is represented by the Flory equation, denoted ex plicitly in the superscript: 2 . 16 Appendix A conditions. First, that the packing factor, whi ch is the fraction of the total system, is universal. Second, the molar volumes are assumed to be additive on reaction , meaning that which is satisfied if mass density of reaction products and reactants are all the same as shown in Appendix B . The calculation results in 2 . 17 where with representing the molar volume of species i . Note is independent of concentration. Thus , we conclud appearing in Eq. ( 2 . 12 ) and that 18 2 . 18 Because is indepen dent of composition, dependence of on composition is coupled to composition dependence of . We expect to be independent of composition when is constant, defined by Flory as an athermal solution ( , and approximately independent of composition when residual contributions are small. Components often have large differences in molar volumes (thus molar densities) creating a composition dependence of molar volume (thus molar density). However, mass den sities, , are frequently found experimentally to be in the range of In cases where the reactants and products of a reaction have the same mass densities, the molar volumes are additive ( ) as shown in Appendix B . , will be constant. Note that the Flory configurational contribution to cancels the molar density appeari ng in Eq. ( 2 . 12 ). Also, this analysis indicates that when residual activity coefficients are used with , the activity coefficient models should not include the combinatorial contribution to avoid double - coun ting the combinatorial effect. 19 APPENDICE S 20 Appendix A : Development of the activity - based equilibrium constant for athermal solutions Here we develop the expression for [30] for the activity coefficients. Considering the following expression: A. 1 where is the volume fraction for component i defined as: . Therefore , the equilibrium constant ratio for products and reactants is A. 2 If the molar volumes are assumed to be additive on product formation and , and recognizing the molar density the e quilibrium constant ratio can be written as: A. 3 A. 4 A. 5 where we define for convenience a composition - independent term A. 6 21 Appendix B : Effect of mass density on volume of mixing H ere we show that if mass density is constant upon mixing or reaction that excess volume is zero. If mass density of the mixture and components are all equal: B. 1 where and is the molecular weight and molar volume of component i respectively. Using the definition of average molecular weight , and selecting E as a basis for the right - hand side: B. 2 Rearrange for B. 3 Similarly, if the mass densities are equal, from Eq. ( B. 1 ) : B. 4 Substituting Eq. ( B. 4 ) in Eq. ( B. 3 ), shows that the excess volume is zero B. 5 Because most liquids have very similar mass densities, the excess volume can often be neglected in all but the most precise work. 22 Chapter 3. Relation of Wertheim A ssociation Constants to Concentration - based Equilibrium Constants for Mixtures with Chain - forming Components Several association modeling approaches have been developed to accurately describe the properties of polar solutions. Chemical theory and Wertheim most popular of these and they have been shown to yield similar functional forms for the contributions of association to Helmholtz energy and activity coefficients. In this paper, we study Flory polymerization theory th rough the work of Campbell and elucidate its correlation to - based equilibrium components h ave one acceptor and one donor site. Algebraic and numerical proofs are given for fluids and mixtures. Additionally, a new generalized activity expression is devel oped for 3.1. Introduction 3.1.1. Motivation As global demand for renewable chemicals and fuels increases, new pathways for their manufacture are being proposed. These projects include upgrading of natural oils, and fermentation of pyrolysis product s to chemicals and fuels. The ubiquitous presence of polar components in biobased process streams complicates the design of industrial separation units. Traditional thermodynamic models, such as non - random two liquid (NRTL) and UNIversal QUAsiChemical (UNI QUAC), are incapable of accurately representing LLE and VLLE in these systems without utilizing more parameters than can be determined with confidence of physical relevance [31] . 23 The limitations of traditional models for polar systems are largely due to their crude representation of important molecular interactions; specifically, hydrogen bonding. Polar molecules that hydrogen bond (or associate) with each other form clusters of chains and/or rings which cause large deviations from ideal behavior. These interactions are responsible for higher boiling points and entha lpies of vaporization of hydrogen bonding compounds compared to molecules of comparable size and polarity that lack the ability to hydrogen bond. The extent and nature of association varies under different conditions including solution composition, tempera ture and pressure. Therefore, it is necessary for thermodynamic models to incorporate hydrogen bonding explicitly rather than lumping the contributions with nonpolar (van der Waals) and simple dipole - dipole interactions. 3.1.2. Historical perspective Associatio n models can be classified into three categories: chemical theories, lattice models and perturbation theories. Apelblat [32,33] presents a historical summary of the advancements made in association mode ling from 1884 to 1984. Among the most noteworthy contributions is [34] theory representing hydrogen bonding as an equilibrium reaction governed by an species formed through hydrogen bonding (e.g. dimers, trimers etc.) to otherwise behave ideally. In other words, physical interactions between the complexes in solution are ignored. This simplification is difficult t o justify and has proved inadequate for many systems especially in modeling liquid - liquid immiscibility. The Dolezalek approach was therefore enhanced by combining the chemical theory contribution with a conventional model, such as NRTL [35] or UNIQUAC [36] , to capture the phy sical interactions. 24 With chemical theory, several decisions must be made a priori regarding the species that exist in solution and the equilibrium constants for various hydrogen bonds. The continuous association model, also known as the Mecke - Kempter mode l, makes the simplification that the value of the equilibrium constant does not depend on the chain or ring size formed. While there has been some discussion on the contradiction of this assumption with spectroscopic evidence [37] , it is often adopted as a firs t - order representation of the complexation in solution. Lattice models represent the system as a lattice of empty and occupied sites and investigate the extent of association by quantifying interactions between adjacent occupied sites. Perhaps the most pro minent lattice model is the one presented by Guggenheim [38] which accounts for non - randomness in the solution. Hydrog en bonding has since been integrated into several lattice models [39 42] . Perturbation t heories add the association to a non - associating reference fluid. Gubbins [43] recently presented a r eview of work done in this area. Several perturbation theories [44 46] have [47 50] , which is the basis of a number of equation of state (EOS) models including Elliott - Suresh - Donohue (ESD) [10,51] , cubic plus association (CPA) [11,52] and the statistical associati ng fluid theory (SAFT) family of equations [9,53 57] . Numerous independent researchers [58 60] have compared association models. Economou and Donohue [58] found that for many systems, both chemical theory (using activity - based and despite differences in their origins and derivations. Wolbach and Sandler derived rela tions between values of the Wertheim association constant for molecules with multiple sites [59,61] . Chapman 25 et al . [9] noti - based ). Campbell rearranged the Flory polymerization theory and was able to find a similar algebrai c form but the terms were not identical [16] polymerization theory has not yet been proven. In this work, we comp are the two approaches for systems of chain - forming species in which the associating components each have an acceptor site and a donor site to contribute to association. Examples of these systems include pure alcohols, alcohol + inerts and alcohol mixture solutions. We also provide algebraic and numerical proofs of universal, and the excess volume is zero. 3.2. Theory and analysis ither deviations from ideal gas mixtures using fugacity coefficients, or deviations from ideal solutions using activity coefficients as discussed by Prausnitz and Tavares [62] . Fugacity coefficients are readily calculated with an equation of state through the compressibility factor, Z . It is common practice to separate contributions to Z relative to the ideal gas value, , according to the nature of the nonideality causing the deviations. 3 . 1 Here and represent effects of the repuls ive and attractive forces and accounts for hydrogen bonding between molecules. Contributions to the Helmholtz energy departure at fixed T and V are correspondingly 3 . 2 In a model such as the SAFT EOS, a bond term is added to both and . An alternative approach is to consider the excess Gibbs energy ( G E ), related directly to the activity 26 coeffici caused by differences in the size and shapes of molecules are accounted for by an entropic combinatorial contribution which ignores all attractive interactions and is the only term remaining at infinitely high temperatures. The association contribution represents non - idealities due to hydrogen bonding and the residual term accounts for other interactions that are not as well - understood and typically includes adjustab le parameters. The contributions combine in a multiplicative form. However, it is conventional to calculate the logarithms of the activity coefficients and present the relationship as 3 . 3 The focus of this work is the term. We begin by introducing the nomenclature used to describe the apparent and true concent overall amount added in composing the mixture divided by the mixture volume. The speciation - dependent. For example, in a pure alcohol, the apparent mole fraction of alcohol is unity, but the solution can also be described in terms of concentrations of true monomers, dimers, trimer s and so on. Notations vary in literature, and the distinction between apparent and true molarity is not always clear. In this work, a subscript T molecular sp ecies which is described as monomer and denoted with the subscript M, which is 3.2.1. Association term in the formation of complexes. Each molecule type in solution may have any number of acceptor 27 ( A ) and/or donor ( D ) sites and bonding between each pair is governed by an association constant, . Notation for the acceptor and donor sites subscr ipts differs in literature, thus comparisons must be done carefully. In early literature, the i and j denote the host component. However, as noted by Michelsen [63] , the summations are gre atly simplified by indexing the acceptor and donor sites rather than the site hosts because site hosts can have more than one type of each site. We write summations over sites in this work. The acceptor and donor terminology allows the model to lend itself to any type of complexation, including covalent bonding. Depending on the number and types of sites attributed to a molecule of a component, different association schemes can be used as reviewed by Kontogeorgis and Folas [64] . Alcohols, for example, are commonly modeled using the 2B association scheme which corresponds to having one acceptor and one donor site per molecul e ( Figure 3 . 1 ). We follow this convention in the present work. A detailed explanation of [9] and Zmpitas and Gross [65] . Figure 3 . 1 : 2B a sso ciation scheme. Hydrogen bonded oligomers with two bonding sites per molecule. The acceptor and donor sites are labeled on the leftmost molecule. The Helmholtz energy departure at fixed temperature and volume due to association derived from eory is implemented in SAFT and CPA as 3 . 4 where is generic notation for a n acceptor, donor, or bivalent (acceptor/donor) site, is the apparent number of moles of the component hosting site , and is the number of identical 28 sites of type B on the host. The term is the fraction of sites of type B that remain non - bonded at equilibrium, calculated through the balance on the probability of binding: where is generic notation for a site which mates with site , denotes the component apparent mole fraction of the site host and is an association constant parameter between site i of type B and site j of type E . In this work is molar density and thus has units of molar volume. With chain - forming components, when each hydrogen bond f ormed, one acceptor ( A ) and one donor ( D ) site are consumed. Thus, in a pure fluid with the 2B association scheme, the fraction of molecules non - bonded at each acceptor site is the same as those free at the donor sites. Therefore, numbering acceptor and do nor sites on the same host using the same subscripts, For mixtures of associating species with the 2B model, Eq. ( 3 . 6 ) is contingent on , which is true when certain cr oss - coefficients equations are used. The reader should recognize that site indexes i and j match the host indexes for the 2B model indexing discussed here. The combining rules we use to calculate cross - coefficients are consistent with Eq. ( 3 . 6 ) and are discussed in Section 3.2.6 . To derive a formula for Hendriks [66] in deriving a generalized function which reduces to Eq. ( 3 . 4 ) when maximized an d results in the association contribution to chemical potential. This approach is extended to calculate the fugacity coefficient of each component in the mixture in Elliott and Lira [8] (Important 3 . 5 3 . 6 29 corrections are made in the errata [67] ). Calculating the activity coefficient from the fugacity coefficient for a given component for the case of conventional mixing rules, van der Waals , a universal packing factor, and no excess volume yields ( Appendix C ) where is given by Eq. ( 1 . 4 ), and represents the fraction of sites free when the host for the site is pure. These constraints are i activity coefficient equation is easier to apply than the equation of Fu and Sandler [6] because derivatives o f are not needed. A more general expression is given in Appendix C , Eq. ( C. 14 ). 3.2.2. - based equilibrium constant In chemical theory, the energetic effects associated with hydr ogen bonding are embedded within an equilibrium constant, . This parameter could be based on the concentrations ( ), mole fractions ( ) or activities ( - known work on oligomer modeling used a concentration - based equilibrium constant and provided a discussion on the inconsistency of with statistical mechanics [30] and variations of it have since been applied in several publications. Chemical kinetics are commonly described in practice using concentrations, and for reversible reactions are consistent with , which is frequently assumed to be independent of composition. This simplification is explored in Chapter 2 where we prove the assumption to be valid only when is invariant with composition. For pure alcohols or alcohol + inert systems, some of the variations of for the formation of a chain of size i+1 from a complex of size i and a monomer are shown in Table 3 . 1 . Here, represents 3 . 7 30 the pure alcohol molar density. The quantity is the volume fraction. The variable , is the true molar density which is smaller than the apparent molar density , because bonding redu ces the number of moles, . In Eq. ( 3 . 11 ), r alc is a measure of the molecular size of the associating species relative to some reference ( ) , chosen by Campbell to be the molar den sity of methanol at 303.15 K. The term denotes the true mole fraction of component i . Table 3 . 1 : Concentration based equilibrium constants . Variations of the concentration - based equilibrium constant u sed in literature. The subscript numeral on differentiates between models. References Eqn Flory [30] , Renon and Prausnitz [68,69] (eqn. 2) 3 . 8 Flory [30] , Nagata [70] , Nath and Bender [71] , Renon and Prausnitz [68,69] (eqn. 1) 3 . 9 Nagata [72] 3 . 10 Brandani [73,74] , Campbell [16] 3 . 11 3.2.3. Relationship between , To understand the relationship between and , we use equations for a pure component with one acceptor and one donor. These are later extende d to mixtures using combining rules in S ection 3.2.6 . We begin with the thermodynamic definition of the activity - based equilibrium constant: 3 . 12 where is the standard state pressure, typically 1 bar. Inserting Eq. ( 3 . 8 ) into Eq. ( 3 . 12 ) yields 3 . 13 31 Economou and Donohue [58] discuss different conventions used for equilibrium constants. By comparing the monomer fraction derived from various models, they conclude that must be proportional to at the same T . Wolbach and Sandler [59] explore the relationship further and show that in the case where the vapor phase is an ideal gas where is a constant that depends on the association scheme used to model the component in solution. For liquid components modeled with 2B, Eq. ( 3 . 14 ) becomes where is a function that depends on the equation of state used to calculate the fugacity coefficient ratio. The value of this function d epends on the mixing rules used for the true solution of monomers and oligomers. The most accepted and widely used combining rules for complexes, introduced by Heidemann and Prausnitz [75] , are 3 . 16 3 . 17 where and are the interaction and size parameters for monomers respectively. Heidemann and Prausnitz [75] show that when Eqs. ( 3 . 16 ) and ( 3 . 17 ) and the Lorentz - Berthelot mixing rules ( , ) are followed, then g is given by Eq. ( 3 . 18 ). Here, is the packing factor. From this equation, an important observation can be made. While of Eq. ( 3 . 15 ) is also a function of the density of the solution due to the function. Kontogeorgis and Folas [64] summarized the 3 . 14 3 . 15 3 . 18 32 assumptions taken in various models in literature and the resulting forms of . When the van der Waals repulsive term is used, Eq. ( 3 . 18 ) gives the radial distribution function at contact distance, . Moreover, applying the combining rules (Eqs. ( 3 . 16 ) and ( 3 . 17 )) to the van der Waals equation for the component fugacity coefficients provides a useful relationship between the fugacity coefficient and true compressibility factor, , as shown in Appendix D . Calculating the ratio of the fugacity coefficients for the reaction of a complex of size i and a monomer to form a complex of size i+1 yields 3 . 19 or 3 . 20 Substituting Eq. ( 3 . 20 ) into Eq. ( 3 . 15 ), then converting to apparent density and inserting Eq. ( 3 . 13 ): 3 . 21 Recognizing that yields a remarkably simple - relationship for pure flui ds based on the definition of Eq. ( 3 . 8 ) . 3 . 22 This result is intuitive considering the conclusions arrived at by several investigators [58,60,66] powerful in that it allows one to use and compare both association models with ease. 3.2.4. Association term Pradhan et al. [76] developed a chemical theory model by integrating an equilibrium constant into the Flory - Huggins (FH) theory. Campbell [16] extended the expression to systems containing any 33 number of alcohols. In this section, we clarify the relationship between this set of equations and other association models. Following Flory [30] , Campbell writes where and are the dimensionless true monomer concentrations of component i in a mixture a nd at purity, respectively. All concentrations , are made dimensionless by dividing by , the liquid molar density of methanol at 303.15 K true complexes in solution which Campbell showed is given by where the summations are over apparent species and is the apparent dimensionless concentration of component i . As show n in Appendix E can be recognized as the total molar concentration of species with a free end group i . The Flory concentration - based is the equ ilibrium constant for the association between a chain of length n with a non - bonded acceptor i, forming a bond with monomer j donor, resulting in a chain of length n + 1 with a non - bonded j acceptor end. Campbell does not discuss acceptors or donors, but t he nomenclature helps clarify how, when the alcohol monomer donor bonds to an end of a chain, it can still participate through its free acceptor. Note that the nomenclature for acceptors and donors could be flipped without changing the balances. Campbell p rovides mathematical derivations for Eqs. ( 3 . 23 ) and ( 3 . 24 ). The free end on a chain could be an acceptor or a donor. Considering the end to be an acceptor, represents the dimensionless molar conce ntration of non - bonded acceptors of type i . Because each alcohol 3 . 23 3 . 24 34 molecule possesses one acceptor site, the dimensionless apparent molar concentration of all acceptors on component i in solution is equal to . Therefore, the fraction of unbonded acc eptors is or . Substituting the dimensionless concentration in terms of molar densities , we find that Therefore, can be written as . Simple manipulation will show that the right expression of Eq. ( 3 . 24 ) is the same as Eq. ( 1 . 4 ). 3.2.5. Algebraic e quivalence of c t heory For a binary alcoho 3 . 7 )) gives and by combining Eqs. ( 3 . 24 ) and ( 3 . 26 ): 3 . 29 3 . 30 3 . 25 3 . 26 3 . 27 3 . 28 35 Next, we calculate the values of and . Campbell provides the relationship between the concentration of monomers ending with group i and given as . Substituting into this equation with the definition of and Eq. ( 3 . 28 ), we find for both components in the mixture and at purity 3 . 31 3 . 32 Theref ore, by combining Eqs. ( 3 . 29 - 3 . 32 ), Eq. ( 3 . 23 ) for the alcohol becomes 3 . 33 Adding and subtracting 1, and rearranging yields 3 . 34 [30] for the combinatorial contribution to activity coefficients is 3 . 35 Similarly, rearranging Eq. ( 3 . 23 ) and substituting Eqs. ( 3 . 29 - 3 . 32 ) for the inert solvent gives 3 . 36 3 . 37 36 3.2.6. Calculating association parameters To clarify the c omparison further, we provide some numerical calculations. Early work done in determining for binary alcohol + inert systems used data reduction of mixture data to fit association constants [68,69] . Commonly, the value of is assumed to be independent of solvent species or concentration. Therefore, models based only on the properties of the pure associating compound were developed. Nath and Bender [71] based on the temperature and enthalpy of vaporization at the normal boiling point of the alcohol. Lobien [77] calculated sm aller values than others using limiting activity coefficient data. Because these approaches differ considerably, there is significant scatter in the resulting values for the same alcohol [71] . For CPA and SAFT, the values are typically determined by fitting density and vapor pressure [9,11,78] . In this work, to provide numerical calculations, we follow the approach outlined by Nagata [70] in applying the Brandani [73] model essentially assumes that the only differences between the vapor pressure of an alcohol, , and that of a non - associating compound with an identical molecular weight and chemical formula, i.e. its homomorph, , are the consequences of hydrogen bonding. Therefore, is determined by matching the two vapor press ures through Eq. ( 3 . 38 ). 3 . 38 where is the free contribution to the second virial coefficient evaluated using Hayden and [79] and is the activity coefficient of the alcohol monomer in pure The volume and mole fractions of monomer in pure alcohol are calculated as [70] 37 3 . 39 3 . 40 Ethers were selected as homomorphs for alcohols and their vap or pressures were determined using [80] correlation with the same n* (ef fective number of carbon atoms in the ether homomorph) values as Nagata. Our resulting values of are given in Table 3 . 2 . Table 3 . 2 : Equilibrium constant for primary alcohols . and calculated at 50 °C Alcohol Methano l 111.0 4549 Ethanol 81.48 4803 n - Propanol 71.27 5381 n - Butanol 69.88 6458 n - Pentanol 57.09 6233 We have refitted the values to be consistent with all of our physical properties. We assume [70] values for the enthalpy for the formation of a hydrogen bond. To extend the analysis to alcohol mixtures, the cross - association parameters for dissimilar components are related through combining r ules. Common combining rules that have been developed and tested [11,51,54,61,81 83] are variations of an arithmetic or geometric mean of the gy and volume parameters. Kontogeorgis [64] outlines a number of these and provides a brief summary of work done in ide ntifying the most accurate rules perturbation theory yield identical results for alcohols described with a 2B model, regardless of the combining rules used if Eq. ( 3 . 22 ) is maintained for the cross coefficients. Therefore, further exploration of the effect of the combining rules lies beyond the scope and intentions of this work and for alcohol mixtures, we use Suresh and Ellio [51] combining rule given by: 38 3 . 41 3.2.7. We implement the abovementioned approaches to provide numerical ev idence for the equivalence approaches produce exactly the same results for systems in which the associating components can be described by the 2B association s cheme. This is illustrated in Table 3 . 3 for two alcohol + inert systems at 30 and 80 mol% of alcohol and one alcohol + alcohol system at 30 mol% of methanol. Both theories yield the same monomer fractions and activ ity coefficients. the perturbation theory and the correct re lationship between the association parameters of the two Future wo rk in this area will look at extending the analysis to components with multiple acceptors and donors. 39 Table 3 . 3 ems. . The reference 3 . For methanol + ethanol, the cross coefficient is 2559 cm 3 /mol or = 62.44 System Ethanol (1) + n - heptane (2) n - Pentanol (1) + n - heptane (2) Methanol (1) + Ethanol (2) 0.3 0.8 0.3 0.8 0.3 T (K) 344.9 344.7 368.4 383.1 347.0 Component i 1 1 1 1 1 2 Chemical Theory 67.37 67.65 53.60 40.30 60.76 64.16 0.2680 0.1790 0.2427 0.1419 0.1025 0.0944 0.0942 0.0749 0.0849 0.1203 0.0966 0.0101 0.0123 0.0119 0.0184 0.0043 0.0095 0.0128 0.012 8 0.0149 0.0192 0.0145 0.0134 0.0301 0.0230 0.0396 0.0527 0.0056 0.0124 0.7195 0.06100 0.5280 0.03109 - 0.002853 - 0.0004456 Wertheim's Theory (cm 3 /mol) 2761 2773 2197 1652 2490 2629 0.3167 0.1696 0.3633 0.2566 0.1361 0.1330 0.0301 0.0230 0.0396 0.0527 0.0056 0.0124 0.9281 0.09320 0.5513 0.03344 0.02999 0.004305 - 0.2086 - 0.03220 - 0.02 337 - 0.002351 - 0.03285 - 0.004751 0.7195 0.06100 0.5280 0.03109 - 0.0028533 - 0.0004455 40 3.3. Conclusions model and chemical theory and compar e the two models. The focus is for the 2B scheme where each associating component has one acceptor and one donor site per molecule such as alcohol + inert and alcohol + alcohol mixtures. In these cases, the association strength in Werthe theory is found to be equal to the concentration - based equilibrium constant polymerization theory. This elucidates the relationship between the two association models. ical theory and show that the similarities drawn by Campbell between identical. Finally, we provide in Eqs. ( 3 . 7 ) and ( C. 14 ) a new generalized expression for activity coefficients from Wertheim theory useful when the packing factor is universal and excess volume is zero. 41 APPENDICE S 42 Appendix C : Derivation of the Wertheim association contribution to the activity coefficient Following Michelsen and Hendricks [66] and to abbre viate notation, we define, at equilibrium, C. 1 Elliott and Lira [8] show that C. 2 C. 3 assoc . By definition, Helmholtz energy is a natural function of temperature and volume. C. 4 For an activity coefficient, we need a derivative of Helmholtz energy at constant T and P [6] . At constant T, C. 5 The expansion rule gives C. 6 The association contribution to the fugacity coefficient can be calculated from C. 7 and C. 8 And making the appropriate substitutions into Eq. ( C. 6 ): 43 C. 9 Substituting Eqs. ( C. 1 and C. 2 ) yields the universal equation C. 10 To sim plify this further, we assume that is given by the van der Waals equation. Applying Eq. ( 3 . 18 ) results in . C. 11 Eq. ( C. 10 ) becomes C. 12 If excess volume is zero and the packing factor is universal, C. 13 Using this equ ation, assoc can be calculated as C. 14 This results in Eq. ( 3 . 7 ) in the text. Note that when excess volume is zero. 44 Appendix D : Relationship between the fugacity coefficient and true compressibility factor Beginning with the fugacity coefficient relation, w e substitute van der Waals equation of state D. 1 and its implementation into the Helmholtz departure equation D. 2 For a true associated species in a pure fluid D. 3 Applying Eq. ( 3 . 16 ) and ( 3 . 17 ), calculating Eq. ( D. 3 ) for a monomer and a dimer gives D. 4 or D. 5 45 Appendix E : the molar concentration of spe cies Here, we provide proof that is the true concentration of all chains ending with i . Beginning with the definition of [16] paper: E. 1 where E. 2 Multiplying through Eq. ( E. 1 ), E. 3 s of length n ending in i is given by E. 4 This equation provides a recursion for simplification. Recognizing Eq. ( E. 4 ) in the second, third, and fourth term on the right hand side of Eq. ( E. 3 ) can be expressed as E. 5 Then E. 6 Simplifying the third and fourth term on the right with Eq. ( E. 4 ) agai n, E. 7 46 This procedure applied continuously to all sums in Eq. ( E. 3 ) yields the result that is also the true con centration of all chains that end in i. E. 8 47 Chapter 4. Relation of Wertheim and Chemical Theories in Systems with Association and Solvation Hydrogen bonding is usua theory. Although the basis of the two methods is different, several researchers have noted parallels between them. In fact, our previous study derives analytical and numerical proofs of th e work, we extend this finding to a complex ternary system containing a chain - forming component and two non - chain - forming components. The non - chain - forming m olecules have only acceptor sites or donor sites, but not both. One molecule has two acceptor sites, and the other a donor site. The main objectives of this work are to elucidate the complicated counting that must be done in chemical theory to correctly ac count for all the bonded species that can form in solution and to demonstrate the equivalence of the two approaches for this complex system. The work also 4.1. Intr oduction As mixtures containing water, alcohols, organic acids and esters occur frequently in the chemical industry, new methods are being developed to accurately model their thermodynamic behavior. The challenge with these types of components is their pro pensity for hydrogen bonding which has a significant effect on many of the mixture physical properties. For example, hydrogen bonding is responsible for the high melting and boiling points of water relative to other chemical compounds of similar molecular size. To robustly capture this behavior, it is now widely accepted that hydrogen bonding, or association, must be represented explicitly as a separate term in thermodynamic models. 48 Association can be modeled several ways and the most common methods fall i nto one of two categories: chemical theory or perturbation theory. Chemical theory is based on the representation of hydrogen bonding as a chemical reaction, described by an equilibrium constant , which results in the formation of a new species. This concept was originally proposed by Dolezalek [34] and has since been implemented to understand on a wide variety of associating systems including l iquid metals [16,37,68,75,84 91] . In particular, Nagata et al. have published an extensive body of work [84,92 94] applying chemical theory to a wide range of organic mixtures. With alcohols, a simplified model, called the Mecke - Kempter model, is often used which assumes that the value of the equilibrium constant is the same regardless of the size of the cluster formed. A challenge in chemical theory that has yet to be resolved is that the types of species that may be formed in solution must be known, or assumed, a priori and the equatio ns for multicomponent mixtures with arbitrary sites are complex. thermodynamic modeling. Developed in the mid - 1980s [47 50,65] , theory is a statistical mechanics approach derived for molecules with a repulsive core and attractive can be used for any type of association including charge transfer complexes. For hydrogen bonding, acceptor sites are allo cated to the electronegative atom(s) which will accept a proton and donor sites are allocated to the protons. Huang and Radosz [55] tabulate schemes that have been develo ped for site allocations and the common 2B scheme used for alcohols is given in Chapter 3. and repeated here ( Figure 4 . 1 49 Figure 4 . 1 : 2B association scheme. The 2B scheme is commonly used to model alcohols by calculating the association between an acceptor site on the oxygen and a donor site on the hydrogen. The probability of association is described by a n association constant, between each binary pair. [11,52] and the SAFT family of equations [9,54 56,95] . In our previous work [96] , we created an activity Michelsen and Hendricks [66] . We also show that with reasonable simplifications, the association constant is equivalent to the concentration - based equilibrium constant and that the two association theories yield the same results for the activity coefficient for alcohol + inert and alcohol + alcohol systems. This aligns with prior work beginning with Elliott et al. [10] who were the first to document the relationship between the activity - based chemical association constant , K a, and perturbation theories. Economou and Donohue [58] then used this equivalence analytically for pure and binary mixtures in which associating components possess 0, 1 or 2 sites each. Similarly, Campbell et al. [16,85] recognized the same parallels using the concentration - based equilibrium constant algebraically for mixtures containing any number of alcohols and, later, for those containing any number of alkanes and any number of alcoh ols. Perhaps the most significant tracking the cluster size and end molecule identity of each cluster. While this scheme is innovative, time - consuming and cumbersome for more complex mixtures, such as water + alcohol, where the number of sites per molecule is greater. 50 In this work, we outline the calculations in chemical and Wertheim theories for a hypothetical ternary mixture in which on e of the components possesses two identical acceptor sites and no donor sites. The presence of such a component complicates chemical theory because certain species resulting from hydrogen bonding can be formed in multiple ways. The main objective of the cu with that of chemical theory. To this end, we implement a counting scheme and show that the monomer fractions obtained by the two association theories are identic al. 4.2. Theory and a nalysis The hypothetical system chosen for this work consists of a component with two identical acceptor sites, such as 1,4 - dioxane, a component with one donor site, such as chloroform (where the electronegativity of the chlorine atoms make s the hydrogen atom susceptible to hydrogen used to model alcohols. We use E, F, G to designate components rather than A, B, C to avoid confusion in use of A for acceptor and C for concentration. The components are illustrated in Figure 4 . 2 . Figure 4 . 2 : Hypothetical molecules for analysis. The system studied in this work i s a ternary mixture of these three molecule types, where the acceptor sites are solid dots and the donors are open dots 4.2.1. p erturbation t heory At equilibrium, molecules which hydrogen bond form clusters of various geometries and sizes. 51 molecules that remain as monomers at equilibrium is the key ter m on which all subsequent calculations depend. Therefore, it is of particular importance that this quantity is known with certainty. When calculating mole balances, a distinction is important between a balance which counts individual molecules as if no ass ociation is occurring, the apparent mole balance, and that which counts each species formed after association has taken place in solution, the true mole balance. The number of species decreases when bonding occurs. Therefore, the true total mole count, , is always less than the apparent mole count, . the fraction of free sites, , is used. This represents the fraction of sites of type that remain non - bonded at equilibrium, where is the number of moles of site that are free and is the total number of moles of site in solution . Variable i can be an index or can be a site host identity; in this work we use the site host identity because the two acceptor sites on E are equivalent. B is a generic variable that can represent an acceptor A , donor D , or carboxylic acid site C which de scribe a site that can bond with other C to form dimers. is calculated via where is generic site notation that can represent an acceptor, donor, or C - site to which can bond. Variable is the apparent mole fraction of the component which hosts , is the number of sit es identical to on the host for . molar density and is the association constant for an interaction between and . 4 . 1 52 Applying Eq. ( 1 . 4 ) to the four unique sites in the hypothetical system yields the equations in Table 4 . 1 . Table 4 . 1 : Equations for sites in hypothetical system . The eq uations are used to calculate the fraction of free sites of each type Interactions involving Wertheim Equation Eq A E, Acceptors on E 4 . 2 A G, Acceptor on G 4 . 3 D G, Donor on G 4 . 4 D F, Donor on F 4 . 5 Here, , , and are the free site fractions of the acceptors on component E and G and donors on G and F respectively. For a given molecule of component E , the monomer fraction is equal to the joint probabilities that the first acceptor is free, and the second acceptor is free. If interactions at the two sides are considered to be independent of one another and recognizing tha t represents the probability that site is free 4 . 6 wher e is the moles of E monomer, is the total apparent moles of E molecules and the represents probability. Similarly, for components F and G : 4 . 7 4 . 8 4.2.2. Implementation of c hemical t heory Calculating monomer fractions with chemical for two reasons. First, chemical theory is a species - wise approach rather than a site - wise approach. Therefore, one must have a priori knowledge of all the possible species that may form in a solution of components. Figure 4 . 3 shows the bonding schemes that can form in the hypothetical ternary 53 mixture studied here. Because molecules of component G possess both an acceptor and a donor site, they can form chains in solution. Therefore, it should be noted that each molecule of component G in Figure 4 . 3 can be replaced by a chain of G molecules. Figure 4 . 3 : Bonding schemes for h ydrogen bonding in a mixture of E , F and G . The bonding schemes circled and denoted by * can be formed in two ways. by four key types of association upon whic h the remaining associations are related. Each can be described by an equilibrium constant as shown Table 4 . 2 . Table 4 . 2 : Key associations in the hypothetical system a nd their equilibrium constants . 4 . 9 4 . 10 4 . 11 4 . 12 A subscript of T on concentration indicates a true value. The common assumption is that the change in Gibbs energy for each bond type is independent of cluster size, thus logarithms of the equilibrium constants add, and the for any compound is the product of the 54 constituent species. From these four interactions, the hydrogen bonded clusters can be formed as summarized in Table 4 . 3 . Table 4 . 3 : Other species formed in hypothetical system . This table outlines the equilibrium constants associated with each species formation and the relation to the key a ssociation constants. 4 . 13 4 . 14 4 . 15 4 . 16 4 . 17 4 . 18 4 . 19 4 . 20 4 . 21 4 . 22 The second complexity with chemical theory is correctly accounting for species that can be formed in multiple ways and which of the ways are indistinguishable by symmetry. In our hypothetical system, all species formed through an unsymmetrical bonding scheme, denoted by an asterisk in Figure 4 . 3 , can be formed in two ways as illustrated for species EF in Figure 4 . 4 . This occurs because E has two identical acceptor sites. To account for these multiple pathways, the equilibrium Table 4 . 2 and Table 4 . 3 . 55 Figure 4 . 4 : Counting species in chemical theory system . Some species may be formed through multiple pa thways. Next, we outline a careful implementation of chemical theory to a mole balance of the three components in our system. To begin with, the mole balance for component E is given by 4 . 23 The double sums in Eq. ( 4 . 23 ) are simplified to a converging series with a closed solution. This calculation is demonstrated in Appendix F for . Substituting the equilibrium constant relations in Table 4 . 2 and Table 4 . 3 for the species concentrations and putting all val ues in terms of and : 4 . 24 Similarly, the mole balance for F and G can be calculated and rearranged as 56 4 . 25 4 . 26 4 . 27 4 . 28 The simplification of the double sums in Eqs. ( 4 . 25 ) and ( 4 . 27 ) is demonstrated in Appendix G . 4.2.3. Numerical e quivalence of c t heories To investigate parallels between the two association theories, it should be recognized that monomer fractions of each component are related by 4 . 29 4 . 30 4 . 31 where the left and right - and the mixture density is calculated by assuming no excess volume of mixing, . 57 B ased on the simplifications of our earlier work [96] , the concentration - based equilibrium constant, as defined in Table 4 . 2 and Table 4 . 3 , is equivalent to the association strength parameter, . Therefore, 4 . 32 For this study, the chosen vales of the association constants do not represent values expected from experiment. Rather, we have chosen parameter values that yield concentrations which are visible on the plot axes (see Table 4 . 4 ). Standard liquid densities of dioxane, chloroform and methanol were adopted for components E , F , G respectively as summarized in Table 4 . 4 . Table 4 . 4 : Association parameters and component densities . Concentration based equilibrium constants are equivalent to Association Parameters (cm 3 /mol) 5000 10000 7000 12000 Component Densities (mol/cm 3 ) 1.167E 02 1.242E 02 2.479E 02 Next, the values of the monomer fraction of each compone nt is calculated. For both association theories, this involves solving simultaneous equations using successive substitution or other iterative techniques. With chemical theory, Eqs. ( 4 . 24 ), ( 4 . 26 ) and ( 4 . 28 ) are solved simultaneously to find values of , and , the equations in Table 4 . 1 are solved simultaneously to yield values for , , and . Bot h methods yield identical results and this is demonstrated for a solution composition of , = 0.2 and =0.6 in Table 4 . 5 and Table 4 . 6 . The sele cted composition provides an environment where complex formation is extensive and only small fractions of the available sites 58 are free. Because the number of species is infinite, the chains are summed in Table 4 . 6 by bonding scheme. Table 4 . 5 : Calculated monomer fraction values . Component monomer fractions for a solution composition of , = 0.2 and =0.6 where = 1.741E 02 mol/cm 3 Chemical Theory Monomer Fraction Value Monomer fraction of E 7.50E 02 Monomer fraction of F 2.92E 02 Monomer fraction of G 5.55E 03 provide species distributions. Table 4 . 6 demonstrates how species distributions can be obtained from both methods. 59 Table 4 . 6 : , = 0.2 and =0.6 . Species Chemical Theory 10 4 (mol/cm 3 ) Eq. ( 4 . 24 ) 2.61 Eq. ( 4 . 26 ) 1.02 Eq. ( 4 . 28 ) 0.58 1.38 2.65 0.67 1.68 5.04 2.56 6.15 3.12 2.43 5.94 3.62 60 Next, we explore the species distribution in solution by varying the composition of the ternary mixture. For t he plots shown in Figures 4 . 5 - 4 . 8 , is fixed at 0.2 (with the exception of the pure G composition point), is varied from 0 to 1 and the balance of the compos ition is . Figure 4 . 5 : True monomer fraction of components in system of study and fractions of acceptors and donors . The mole fraction of component F is fixed at 0.2. Figure 4 . 6 : Distribution of species with component E when =0.2 . The bottom plot is a zoomed panel of the top plot 61 Figure 4 . 7 : Distribution of species with component F when =0.2 . The bottom plot is a zoomed panel of the top plot Figure 4 . 8 : Distribution of species with component G when =0.2 . The bottom plot is a zoomed panel of the top plot 62 Figure 4 . 5 shows the mole fraction of monomers for each component and the percentage of acceptors and donors in solution. The successive three plots show how components E , F and G are distributed in solution. The fraction of acceptors and donors in solution are calculated by the ratio of acceptors (or donors) to the total sites in solution. Acceptors are in excess to the left of and because the association constants are large, almost all of the donors are bonded until above where some F remains free. In Figure 4 . 6 the true mole fraction of monomer E, ( ), stays near the apparent mole fraction of E on the left side of the figure because, although the true moles of E decrease due to bonding, the true mole fractions are calculated using the total moles of all species after bon ding, which is less than the apparent moles, ( or ), used to calculate the apparent mole fraction ( ). For example, at the composition and , , and . This yields a value of 0.766. Moreover, almost all the F is associated, and the dominant species are E and EF. At the apparent E and apparent G are equimolar, so the compounds are favored. At , the ratios and are equal, and thus the dominant species are , , and . The figures show the sum of mole fractions for oligomers rather than oligomer distributions. For a pure sol ution of a chain - forming compound such as component G , it can be shown that the true mole fraction of monomers is always larger than the mole fraction of any oligomer species. The number of moles of a chain of G molecules of length i is [8,67] : 4 . 33 63 where is the total true number of moles and is the true monomer fraction of G . Note that for pure component G , and therefore, can be written as . Further, it is shown that: 4 . 34 where is the total apparent number of moles. Combining Eqs. ( 4 . 33 ) and ( 4 . 34 ) yields: 4 . 35 Using this equation, the relationship between and can be calculated to be: 4 . 36 Recognizing that is the monomer fraction in solution, we note an interesting trend. As the value of is increased, molecules are more likely to form bonds and therefore decreases. This occurs such that is always a fraction and is less than . In other words, there are always more monomers in solution than dimers, more dimers than trimers and so on. This t Figure 4 . 9 for pure methanol at mol/cm 3 for four values of . As a reference, the cubic - plus - association (CPA) equation of state value for the for pure liquid methanol with itself is approximately 16000 cm 3 /mol at 298.15 K. More interestingly, as approaches infinity, the distribution of various cluster sizes becomes uniform. 64 Figure 4 . 9 : True mole fraction of methanol chains calculated with 4 values of 4.2.4. Activity c oefficients As further demonstration of the equivalence of Wertheim and chemical theory methods, we illustrate t he calculation of activity coefficients. In our previous work [96] , we demonstrated that the association contribution to the activity coefficients, assuming van der Waals correction to xcess volume, results in: where the superscript 0 represents the pure state for component k. In our previous work [96] , the logarithm operator was missing on Equation A14 and it should have appeared as 4 . 38 The activity coefficient from chemical theory is 4 . 37 = 65 Note that the molar densities are equivalent num erically to the molar concentrations used above in the reaction equilibria. The previous work [96] showed that the two approaches are related by . Table 4 . 7 summarizes calculation of the activity coefficients from the two approaches using the concentrations determined above. Figure 4 . 10 shows the activity coefficient contributions for mixtures with increasing compositions of component G when the mole fraction of F is fixed at 0.2. The association parameters and the pure component densities used in all cal culations are given in Table 4 . 4 . Table 4 . 7 : Calculated parameters and intermediate values for activity coefficient calculations . The selected composition is , and and = 1.741E 02 mol/cm 3 , = 3.889E 03 mol/cm 3 Component 0.2 0.2 0.6 Chemical theory 2.62E 04 1.01E 04 5.80E 05 1.17E 02 1.24E 02 9.38E 05 1.17E 02 1.24E 02 1.53E 03 1.52 2.51 6.65E 02 Wertheim's Theory 0.274 - 0.190 - 2.91E 02 2.93E 02 1.00 - 6.15E 02 - 1.00 6.15E 02 1.43 2.45 1.08E 02 9.17E 02 6.40E 02 5.57E 02 1.52 2.51 6.65E 02 4 . 39 66 Figure 4 . 10 : Calculated values of the activity coefficient contributions . For each component, heory ( ), chemical theory ( ) contributions are shown as a function of the apparent mole fraction of G . The m ole fraction of component F is fixed at 0.2. 4.2.5. Perspective on f urther g eneralization It is worth noting that the system cho sen for this study is intentionally rather simple compared to those encountered in the chemical industry. By definition, components which contain only one type of site, such as components E and F , limit the extent of association in solution. When a molecul e of F binds to an acceptor site, it seals one end of the chain from further association. With mixtures of water, alcohols and/or glycols, hydrogen bonding is extensive and can create cyclic [97] . Though chemical theory can capture cyclic species [13,98] , specifying the species to include in solution is challenging. However, even 67 without considering cyclic species, the intricacy of chemical theory is in stark contrast with the onding is calculated using the site balance introduced in Eq. ( 1 . 4 ). Chemical theory requires correct identification of all possible bonding schemes and the number of combinations through which each can form. It is thus this compl exity researchers [10,16,58] . Indeed, in our own previous work [96] , the algebraic and numerical equivalence of the two theories is shown for simple mixtures containing chain - forming species, specifically pure alcohol and alcohol + inert type systems. However, in those cases, the alcohol, modeled in the same way as component G here, results in only chain equation s that are more easily manipulated. The presence of species E creates a complexity that has not been previously demonstrated. To our knowledge, no preceding work has shown the numerical equivalence of the two methods for multicomponent non - chain - forming sp ecies with this rigor. 4.3. Conclusions hypothetical system that contains a chain - forming component and molecules with multiple sites which can only solvate. For chemic al theory, the calculations are outlined carefully to elucidate the oftentimes complicated counting procedure required to account for all associating species. Moreover, the effect of the value of the association parameters on derived monomer fraction value s is studied. We demonstrated how Wertheim site free fractions can be used to calculate species concentrations. It is found that, when implemented carefully, Wertheim and chemical theory yield identical species distributions for the given system, providing that cyclic species are not present. Because the system of study is a mixture of generalized non - specific components, we expect that 68 the work presented here can be extended to all other chemical systems that do not form cyclic species. Although both theor advantage due to the simplicity of its universal equations and is therefore recommended as the superior method for determining monomer fractions. 69 APPENDICES 70 Appendix F : Simplification of double sums without a multiplication factor Here, we demonstrate double sum simplifications such as that of F. 1 Careful inspection of the equilibrium c onstants shows that whenever sum to the same value, all concentrations where G is redistributed on both sides of E are identical, i.e. F. 2 Therefore, the d ouble sum can be written as a single sum: F. 3 Substituting the definition for the equilibrium constant for this reaction from Eq. ( 4 . 17 ): F. 4 F. 5 Recognizing the c onverging series F. 6 results in the following simple term: F. 7 71 Appendix G : Simplification of double sums with a multiplication factor For the double sums involving an coefficient, such as, G. 1 Applying Eqs. ( F. 2 ) and ( 4 . 17 ) again yields two single sums. = G. 2 The first sum is addressed in Appendix A. The second sum is also a converging series. G. 3 Applying Eqs. ( F. 6 ) and ( G. 3 ): G. 4 G. 5 72 Chapter 5. Applications of a Wertheim association activity coefficient model to methanol and ethanol - containing mixtures A new activity coefficient method for associating systems, Wertheim Association Gamma (WAG), is evaluated for methanol and ethanol - containing binary and ternary mixtures. The method ory and a residual contribution. This residual contribution is represented using several existing thermodynamic models including NRTL, Scatchard - and limitations of the mo dels are compared to those of PC - SAFT and CPA. The WAG - NRTL and WAG - Nagata approaches capture the behaviors of these mixed systems well when modeled with four parameters. Meanwhile, WAG with Scatchard - Hildebrand (WAG - SH) out - performs CPA and PC - SAFT in se veral cases, with only two fitted parameters. Lastly, the limitations of existing methods for calculating association parameters are explored and shown to need improvement. Introduction Primary alcohols have many important uses in the chemical industry. Ethanol has been studied extensively as a promising alternative to petroleum - derived chemicals and fuels [99 101] . Methanol is often used as a raw material in the manu facture of important compounds, such as methyl - tert - butyl ether, and to prevent the formation of gas hydrates in the oil and gas industry [102] . In such cases, the thermodynamic behavior of the systems can be difficult to represent with traditional models and more a dvanced methods must be used. For alcohols, water, organic acids and other polar molecules, it is now well - established that models that include hydrogen bonding effects perform better than those that do not. For this reason, association equations of state, such as SAFT [9] (and its derivatives [54 56,95] ), CPA [11] and ESD [10] provide better representation of these types of systems. 73 Though theoretically more rigorous, equations of state require fitting of several pure component parameters, requiring extens ive vapor pressure and density data that are often not available for new compounds. Lack of data results in large uncertainty in parameter values. Further, the dispersion terms and association terms are both attractive and correct attribution of effects to the respective interactions is difficult when fitting limited macroscopic data. In contrast, activity coefficient ( ) models have been favored by industry for many applications (particularly those at low or moderate pressures) because the fitting of pure - component properties is only required in the range to be modeled, and attention can be focused on the mixture nonidealities. Additionally, the representation of the mixture is not biased by challenges in fitting the pure component behavior. Association h as been modeled in activity coefficient models as early as 1969. Harris and Prausnitz [103] combined the van Laar equati [34] chemical theory to represent physical interactions and hydrogen bonding respectively. The Dolezalek approach describes the formation of a hydrogen bond as an equilibrium reaction guided by an equili brium constant. More recently, Karachewski [13] combined this association theory with Scatchard - Hildebrand [19,20] and derived the value of the equilibrium constants of hydrogen bonding from NMR spectroscopy. Asprion et al. [7] demonstrated that use of chemical theory with local composition models can suppress false prediction of VLLE. An alternative me thod for calculating association is [47 50] . This model has become inc reasingly popular in recent years and is described in more detail in Section 5.2.1 . Fu and Sandler [6] were the first to integrate heory for association into UNIQUAC [ 12] , an activity coefficient model, achieving significant improvement for binary systems. The model was also extended to a predictive UNIFAC method [104] with good results, albeit with the association parameters regressed to experimental phase equilibria data without consideration of actua l spectroscopically measured association. 74 Unfortunately, both of Fu's models include a derivative term that make them difficult to generalize for multicomponent systems. In our previous work [96] , we developed for association which will henceforth be referred to as the Wertheim Association Gamma (WAG) work yielding a general equation that can be applied readily in multicomponent systems. More importantly, the generalized association term can be attached to existing activity coefficient models to enhance their ability to model systems with hydrogen bonding. In this work, we explore the capabili ties and limitations of three variations, WAG - NRTL, WAG - Scatchard - Hildebrand (WAG - SH) and WAG - Nagata, for various systems of industrial importance and then compare their performance to two association equations of state, CPA and PC - SAFT. Model theory and p arameters In this section, we outline the models compared herein. In cases with the WAG models, the vapor phase non - [79] model which has long been recognized f or its accurate modeling of the vapor phase. Liquid phase non - idealities are calculated using an activity coefficient model which is a product of three contributions: 5 . 1 or 5 . 2 The terms , and are the combinatorial, residual and association contributions t o the activity coefficient respectively. The combinatorial contribution calculates the entropic effects that arise due to differences in molecular sizes and shapes in solution. In theory, this would be the only term remaining at liquid densities and inf initely high temperature. The residual contribution accounts for molecular interactions between unlike molecules. In practice, depending 75 on parameters used, this term can remain at liquid densities and infinitely high temperature. The association contribut ion calculates the effects of hydrogen bonding. 5.2.1. Association contribution [47 50] , a statistical mechanics approach that quantifies that probability of bonds forming between sites on molecules. Acceptor and donor sites are allocated on electronegative atoms that accept hydrogen bonds and the protons that can be donated respectively. In our previous work [96] , we introduced a The model assumes conven tional mixing rules, no excess volume, and a universal packing fraction. Further details regarding the development of the model are given in Appendix C of Chapter 3 and are omitted here. The resulting functional form is given by: where and represent the fraction of sites free in the mixture and when the host for the site is pure respectively. is the number of identical sites on a molecule of host component. can be c alculated using: where represents any site that can bond to and is the association constant for the two sites. The variable is the apparent mole fraction of the component that hosts the site. Donor sites are solved by using a similar equation obtained by flipping the acceptor and donor identities in t he equation. 5 . 3 5 . 4 76 5.2.2. Combinatorial contribution The combinatorial contribution is commonly calculated using the well - equation: where is the volume fraction or modified volume fraction for component i calculated by assuming zero excess volume: s original theory (and unless otherwise specified in this work), = 1. Despite its extensive factor. This is reasonable for smaller molecules of similar size wh ere the entropy effects are small. However, the assumption of constant packing factor is inaccurate for polymers in solvents where entropy of mixing is significant [105] . Various attempts to correct such weaknesses in the theory commonly use the Staverman - Guggenheim (SG) theory [38,106] whi ch has been integrated into well - known models including UNIQUAC [12] and UNIFAC [29] . Recent work by Vahid et al. [105] compares the excess entropy term ca lculated using several methods with results from discontinuous molecular dynamics simulations. They found that Flory, both alone or corrected with SG, underestimated the excess entropy for certain values of packing factor and misrepresented its skewness wi th respect to equimolar composition. Other semi - ( 5 . 6 ) with the UNIQUAC parameter and/or using fractions for the power, . Donohue and Pr ausnitz [107] were the f irst to suggest using with a non - unity power in the Flory - Huggins - Hildebrand. Kikic et al. [108] compared the excess Gibbs energy calculation from Donohue and Prausnitz model with others that corrected 5 . 5 5 . 6 77 Flory - Huggins with a modified SG term. They tested both and for the SG term and showed that both modifications far out - p erformed traditional SG for aliphatic hydrocarbons. Following these studies, two groups published modified UNIFAC models that incorporate a non - unity power in the Flory portion of the model. The first, which was developed at the University of Dortmund and is therefore often referred to as the Dortmund modified UNIFAC [109,110] , sets and retains the SG correction. The second, developed at the University of Denmark (Lyngby), is sometimes referred to as the Lyngby modified UNIFAC model [111,112] and follows Kikic et al. . The Lyngby modified UNIFAC model also removed the SG correction completely, citing work by Sayegh and Vera [113] showing that the correction often has little effect or so metimes overestimates the combinatorial contribution leading to unrealistic excess entropies. In our own work, we performed preliminary comparisons between: 1. 5 . 6 ) with =1) 2. 5 . 6 ) with =1) + SG 3. 5 . 6 ) with =2/3) 4. 5 . 6 ) with =1) + van der Waals ent ropy correction [8,67] A survey of several systems showed no appreciable differences in the fitting of binary systems studied in this work. In all cases, the residual term compensated fo r differences in the combinatorial term. Accepting that the original Flory and the Staverman - Guggenheim have some deficiencies, we follow Kikic et al. combinatorial contribution using Eq. ( 5 . 6 ) with =2/3 and no additional correction. The decision to use molar volumes rather than the UNIQUAC parameter was made based on the note by Kikic 78 et al. that the former case showed a slight advantage over the latter in modeling linear alkanes. The use o f volumes provides some temperature dependence in the Flory term. 5.2.3. Residual contribution Three models for the residual contribution are tested in this work: non - random two liquid (NRTL) theory, a model developed by Nagata, and Scatchard - Hildebrand yielding the WAG - NRTL, WAG - Nagata and WAG - SH models respectively. In this section, we provide a brief overview of each of these models. Activity coefficient models for liquid phase modeling may assume random mixing or consider variations in local compositions. In contrast to random mixing models , which assume that species are randomly dispersed in solution, local co mposition models represent the compositional differences between the solvent shell of a molecule and the bulk fluid. This variation is modeled by defining a constant, : where is the mole fraction of molecules of component i around a molecule of component j . represents a ratio of the local composition to the bulk composition of the fluid and is calculated differently in every mode l. The NRTL equation, which was developed by Renon and Prausnitz [3] in 1968, calculates with adjustable parameters, and : Here, is the interaction energy parameter between a molecule of component i and one of compo nent j and is called the non - randomness parameter. The value of is then calculated through: 5 . 7 5 . 8 79 Often, is fitted by two parameters to increase the temperature dependence of the term: NRTL is a commonly used model in the chemical industry because its functional form makes it remarkably flexible . However, its parameters lack the theoretical significance necessary to derive physical meaning or predict their values without regression. In this work, NRTL models (both with and without WAG) are called 2 - parameter model s if is set to 0 and only indicates that there will be two values used to fit each binary system: and . The 4 - parameter model fits both and . The non - randomness param eter, , is not fitted to each individual system. Instead, was set to 0.2 when the system has a liquid - liquid immiscibility and 0.3 otherwise. Another local composition model, Wilson [18] , is the basis of what is referred to herein as the Nagata model. Wilson used t he variable to represent local composition. Ellliott and Lira [8,67] is related to in Eq. ( 5 . 7 ) through: where is a fitted parameter. With careful evaluation of the te [8,67] . For convenience, one can consider Flory to be the combinatorial contribution to the excess Gibbs ener gy, , and attribute the other terms to a residual contribution, . 5 . 9 5 . 10 5 . 11 80 In their 1997 work, Nagata et al. [15] used this idea to combine Eq. ( 5 . 14 ), which they r efer to as molar volume in Eq. ( 5 . 6 ) is replaced by a UNIQUAC volume parameter raised to 2/3, and the chemical cont ribution is calculated in terms of the true monomer activities. This approach is very similar to that developed by Campbell [16] and we show in our previous work that the resulting , is related to the traditional Flory activity coefficient, , through: Theref ore: equilibri where is the volume fraction for component i . The quantity is calculated through Eq. ( 5 . 10 ) and the 2 - and 4 - par ameter implementations of this model are fitted in the same way as described above for NRTL. To avoid confusion, we refer to this as WAG - Nagata, rather 5 . 12 5 . 13 5 . 14 5 . 15 5 . 16 5 . 17 81 than Wilson, to make it clear that only the residual portion of Wilson is used and the association and c ombinatorial contributions are calculated as described in S ections 5.2.1 and 5.2 .2 . The third residual model evaluated here is Scatchard - Hildebrand (SH). Developed in the early 1930s [19,20] , SH is a regular solution theory which extends the van der Waals equation of state to mixtures. The excess Gibbs energy an d have the form: where is the mixture molar volume calculated by assuming no excess volume, , is the volume fraction of component i and The solubility parameter, , is a measure of the internal energy departure per molar volume of a component and is therefore often referred to as the cohesive energy density. of Eq. ( 5 . 20 ) is a binary interaction parameter that can be fitted to experimental mixture data. Because adjusting is equivalent to fitting we regress as: where and are t he fitted parameters. In summary, will be calculated using the following 3 WAG models: 5 . 18 5 . 19 5 . 20 5 . 21 82 5.2.4. EOS model parameters The WAG models developed in Eqs. ( 5 . 22 - 5 . 24 ) are compared to CPA and PC - SAFT. These EOS models are used as a reference and thus, their derivations are not provided here but can be found in references [11] and [95] for CPA and PC - SAFT respectively. We use the implemen tations in version 9 of AspenPlus® and we fit the binary interaction parameter, with the forms implemented in AspenPlus® : 5.2.5. Association parameters Ofte n, of Eqn. 5 . 4 is calculated by fitting to pure component properties. In this work, the WAG models use the CPA association parameter, given by: 5 . 26 where is the CPA covolume term, is an interaction energy parameter, is is a fitted parameter. The term is the radial distribution function at contact given by: 5 . 27 5 . 22 5 . 23 5 . 24 5 . 25 83 The values of CPA parameters used in this work are given in Table 5 . 1 . The values of and in the table refer to the self - association interactions, which occur between acceptors and [55] notation, alcohols are modeled with a 2B scheme (one acceptor site, one donor site) and water with 4C (two identical acceptor sites, two identical donor sites). For ethano l, one set of CPA parameters are used for all the WAG models and another is used for the CPA fits they are compared to. The reasoning for this choice is given in Se ction 5.4 . Table 5 . 1 : CPA parameters used to calculate . The terms and pertain to self - association interactions Compound (K) (L/mol) (bar L 2 mol - 2 ) (bar L mol - 1 ) Reference Methanol 512.6 0.03098 4.053 0.4310 245.91 0.0161 [1] Ethanol (CPA) 513.9 0.04910 8.672 0.7369 215.32 0.0080 [1] Ethanol (Other) 513.9 0.04690 6.701 0.7987 247.72 0.0127 [2] Water 647.3 0.01452 1.228 0.6736 166.55 0.0692 [114] n - Pentane 469.7 0.09101 18.20 0.7986 [114] n - Heptane 540.2 0.1254 29.18 0.9137 [114] n - Decane 617.7 0.1787 47.39 1.1324 [114] Cyclohexane 553.6 0.09038 21.26 0.7427 [114] Cross - association, which occurs between molecules of different components, is approximated by combining self - associati on parameters according to simplifying rules called combining rules. 5 . 28 For WAG models, we implemen t a modification. The ECR includes the radial distribution function of the pure components i and j in and calculates their geometric mean values of the association parameters. However, this eliminates any compositional dependence of for the mixture . Therefore, we modify the rule as follows: 84 5 . 29 Where the pure component values are determined at the pure component densities at the same temperature. Assuming that , Eq. ( 5 . 29 ) can be applied by setting: 5 . 30 5 . 31 where is a constant in this work. For the CPA model, the CR - 1 combining rules are applied for cross - association such that: 5 . 32 5 . 33 For the PC - SAFT model, parameters are obtained from literature for associating [115] and inert [95] components. The association parameter is given by: 5 . 34 where is the radial distribution function for hard spheres calculated with the Carnahan - Starling [116] equation and: where has units of volume and is the tempera ture - independent segment diameter for component i . Eqs. ( 5 . 33 ) and ( 5 . 35 ) are the combining rules used to model cross - association in PC - SAFT. 5 . 35 85 Model comparisons To isolate the capabilities of the residual model for all WAG models, the Wertheim association parameter was calculated using the CPA parameters from Table 5 . 1 consistently for all models except PC - SAFT. Only the residual term parameters were adjusted. Liqui d molar volumes were calculated using DIPPR molar volumes as provided by the Aspen Plus® databases given in Table 5 . 2 . PC - SAFT parameters are given in Table 5 . 3 . Table 5 . 2 : Aspen Plus® databases used in calculations . WAG - * denotes all WAG models Model Databases in AspenPlus ® WAG - * APV90 PURE35, NISTV90 NIST - TRC CPA APEOSV90 AP - EOS, APV90 PURE35, NISTV90 NIST - TRC PC - SAFT APV90 PURE35, APV90 PC - SAFT Table 5 . 3 : PC - SAFT pure component parameters . The terms and pertain to self - association interactions Compound (Å) /k (K) (K) Reference Methanol 1.526 3.230 188.9 2900 0.03518 [115] Ethanol 2.383 3.177 198.2 2653 0.03238 [115] Water 1.066 3.001 366.5 2501 0.03487 [115] n - Pentane 2.690 3.773 231.2 [95] n - Heptane 3.483 3.805 238.4 [95] n - Decane 4.663 3.838 243.9 [95] Cyclohexane 2.530 3.850 278.1 [95] Hayden - by Aspen Plus® ver. 9 (given in Appendix H) without adjustment. Excess volume was assumed to be zero, and all mixture molar volumes used for volume fractions were temperature - dependent. To facilitate regressions, the WAG and residual models were programmed as a user activity coefficient subroutine for AspenPlus®. Data sources used for regression are summarized in Table 5 . 4 and the values of the regressed parameters are given in the Appendix. When multiple data sets are available from the source, all sets are used in the regressions. 86 Table 5 . 4 : Data used for parameter regression . System Data fitted to Methanol + n - heptane VLE [117] , LLE [118] Methanol + cyclohexane VLE [117] , LLE [5] , H E [119] Methanol + n - pentane VLE [120] Methanol + water VLE [121,122] Methanol + ethanol VLE [123 125] , H E [126] n - Heptane + cyclohexane VLE [127] (point at omitted), H E [128] Ethanol + n - decane VLE [129] , H E [130] Ethanol + cyclohexane VLE [131] , H E [132] The goodness of fit for each regression is given in this work by the sum of squares value, SSQ, calculated by: Here, is generic notation for all properties used in the regression such as mole fractions or heats of mixing. Th e superscript refers to the experimental value and is a standard deviation assumed for the experimental data; is 0.1 K for temperature values in the VLE data sets and 0.01 K for those in LLE data sets, 0.1% for excess enthalpy data, 1% for vapor mole fractions, and 0.1% for pressure and liquid mole fraction data. These weights reflect the confidence with which each property can be obtained experimentally. Table 5 . 5 g ives the total SSQ of each model for all systems tested here with the exception of two. The methanol + pentane system was omitted because 4 - parameter WAG models and NRTL were deemed unnecessary due to the satisfactory fits of the 2 - parameter versions. The methanol + cyclohexane system was also omitted because the 2 - parameter WAG - NRTL model did not converge. From these results, the four parameter versions of WAG - NRTL and WAG - Nagata provide the best overall fits. This is not surprising when compared to CPA an d PC - SAFT since the equations of state were fitted with only two parameters each. However, it is encouraging that 5 . 36 87 WAG - SH and the 2 - parameter variation of WAG - Nagata, which are also fitted with just two parameters, perform better than the CPA and PC - SAFT fo r the systems tested. Table 5 . 5 : Total SSQ for systems studied in this work . Systems for which only some model were fitted have been omitted from this table and (2) and (4) denote the number of fitted parame ters Model Number of fitted parameters Total SSQ 10 - 5 WAG - NRTL(4) 4 1.12 WAG - Nagata(4) 4 1.82 NRTL(4) 4 2.13 WAG - SH 2 8.14 WAG - Nagata(2) 2 9.30 CPA 2 11.9 NRTL(2) 2 19.7 WAG - NRTL(2) 2 27.2 PC - SAFT 2 42.5 While the total SSQ provides a reasona ble overview for the fit, it is important to consider each regression individually to understand the limitations of each model. For example, the high SSQ value for PC - SAFT misrepresents the majority of fits. This SSQ is large due to the moderately poor rep resentation of the methanol + n - heptane liquid - liquid equilibrium. To explore fits more carefully, Table 5 . 6 arranges the models according to SSQ value for each system. Each of the models will be represented by the line and/or color in Table 5 . 7 . In some figures, certain lines may be obscured due to overlap with others. 88 Table 5 . 6 : Model performance in order of best fit . Numb ers are SSQ 10 - 4 values and (2) and (4) indicate the number of fitted parameters in the model. (*) denotes that only VLE (not excess enthalpy) data was fitted Methanol + n - heptane WAG - NRTL(4) 8.83 NRTL(4) 13.0 WAG - Nagata(4) 16.0 WAG - SH 75.2 WAG - Nagata( 2) 87.0 CPA 103 NRTL(2) 183 WAG - NRTL(2) 265 PC - SAFT 419 Methanol + cyclohexane NRTL(4) 17.6 CPA 20.6 WAG - Nagata(4) 23.8 WAG - NRTL(4) 24.5 WAG - SH 41.8 PC - SAFT 80.3 WAG - Nagata(2) 107 NRTL(2) 285 Methanol + n - pentane WAG - SH 0.170 WAG - Nagata(2) 0.215 WAG - NRTL(2) 0.241 PC - SAFT 0.279 CPA 0.303 NRTL(2) 0.597 Methanol + water WAG - Nagata(4) 0.133 NRTL(4) 0.169 WAG - SH 0.199 NRTL(2) 0.218 WAG - Nagata(2) 0.244 CPA 0.325 PC - SAFT 0.525 WAG - NRTL(4) 0.564 WAG - NRTL(2) 0.790 Methanol + ethanol WAG - NRTL(4) 0.00872 NRTL(4) 0.00882 WAG - Nagata(4) 0.0108 WAG - NRTL(2) 0.0423 NRTL(2) 0.0431 WAG - Nagata(2) 0.0716 PC - SAFT 0.0993 CPA 0.618 WAG - SH 0.664 n - Heptane + cyclohexane WAG - Nagata(2)* 0.0152 NRTL(2)* 0.0152 WAG - NRTL(2)* 0.0153 WAG - Nagata(4) 0.0172 NRTL(4) 0.0191 WAG - NR TL(4) 0.0193 WAG - SH 0.0881 PC - SAFT 0.0893 CPA 0.101 Ethanol + n - decane WAG - NRTL(4) 1.51 WAG - Nagata(4) 1.68 PC - SAFT 4.35 WAG - SH 4.89 WAG - Nagata(2) 5.23 WAG - NRTL(2) 5.64 NRTL(4) 7.50 NRTL(2) 12.0 CPA 13.0 Ethanol + cyclohexane WAG - NRTL(4) 0.305 WAG - Nagat a(4) 0.318 WAG - SH 0.392 WAG - Nagata(2) 0.401 PC - SAFT 0.443 NRTL(4) 0.560 WAG - NRTL(2) 0.645 NRTL(2) 1.43 CPA 1.98 89 Table 5 . 7 : Legend for plots in Model Capabilities section Model Notation in Figures WAG - NRTL (2 parameters) or WAG - NRTL (4 parameters) or WAG - Nagata (2 parameters) or WAG - Nagata (4 parameters) or WAG - SH or NRTL (2 parameters) or NRTL (4 parameters) or PC - SAFT or CPA or In the upcoming sections, we show phase diagrams, hea t of mixing plots (when applicable), activity coefficients and percent errors for the K - ratios calculated by the various models. As an additional measure of the model fits, the percent error of the K - ratios is calculated for each binary as: where KVL is the K - ratio in the vapor - liquid equilibrium (VLE) region and the superscript of indicates experimental values. For the liquid - liquid eq uilibrium (LLE) region, KVL is replaced by the K - ratio in the LLE region. 5.3.1. Binary system with no association First, we test the capabilities of the residual and combinatorial activity coefficients with a simple binary mixture of two inert components, n - hep tane and cyclohexane. The results are given in Figures 5 . 1 - 5 . 4 . The 2 - parameter WAG and NRTL models were unable to simultaneously fit the VLE and heat of mixing data well. Because the heat of mixing is small and inconsequential for nearly ideal systems such as this, the parameters were fitted only to the VLE data. Therefore, the lines for the 2 - parameter WAG and NRTL models are shown in the VLE diagram but omitted from the h eat of mixing diagram. All the models were able to describe VLE behavior equally well. 5 . 37 90 Figure 5 . 1 : Phase equilibria diagram for n - heptane + cyclohexane at 101.3 kPa . Figure 5 . 2 : Heat of mixing for n - heptane + cyclohexane at T=298.14 K and 101.3 kPa . The WAG - NRTL(4) line is obscured by the WAG - Nagata(4) and NRTL(4). The WAG - SH, CPA and PC - SAFT lines are overlapped. 91 Figure 5 . 3 : Calculated and experimental activity coefficients for n - heptane + cyclohexane at 101.3 kPa . Data and fits correspond to the experimental data in Figure 5 . 1 Figure 5 . 4 : Percent errors of K - ratios for n - heptane + cyclohexane system . Figures (a) is the K - ratio of n - heptane and (b) is the K - ratio of cyclohexane in the VLE region at 101.3 kPa (a) 92 5.3.2. Binary systems with self - association Figures 5 . 5 - 5 . 8 show fits for the methanol + n - heptane system. The fits for heat of mixing exhibit a linear segment that spans the predicted LLE for each model. NRTL and WAG - NRTL(2) ar e fitted phase equilibria and the predicted heat of mixing in Figure 5 . 6 , can be explained by its lack of association . Including the Wertheim association term in WAG - NRTL improves the fit significantly. However, WAG - NRTL(2) does not have the temperature dependence necessary to fit the LLE region accurately and two more parameters must be included. PC - SAFT also shows a po or fit in the liquid - liquid region. The other models are able to represent the phase equilibria with comparable accuracy. Figure 5 . 5 : Phase equilibria diagram for methanol + n - heptane at 101.3 kPa . Experi mental VLE and LLE data is shown by markers. 93 Figure 5 . 6 : Heat of mixing for methanol + n - heptane at T=298.15 K and 101.3 kPa . Experimental data [133] is shown for the miscible region ( < 0.1464 and > 0.8899). Figure 5 . 7 : Calculated and experimental [4] activity coefficients for methanol + n - heptane at 101.3 kPa . Data and fits correspond to the experimental VLE data in Figure 5 . 5 94 Figure 5 . 8 : Percent errors of K - ratios for methanol + n - heptane system . Figures (a) and (c) are K - ratios in the VLE region and (b) and (d) are K - ratios in the L LE region The fitting results of the methanol + cyclohexane system are very similar to those of methanol + n - heptane, which are given in Figures 5 . 9 - 5 . 12 . NRTL and the tw o - parameter WAG models provide poorer fits, particularly in the LLE. The excess enthalpy data in Figure 5 . 10 is fitted well by the four parameter WAG models. While capable of capturing the magnitude well, CPA and W AG - SH are slightly inferior in capturing the minor asymmetry in the heat of mixing. (a) (b) (c) (d) 95 Figure 5 . 9 : Phase equilibria diagram for methanol + cyclohexane at 101.3 kPa . Experimental VLE and LLE data is shown by markers. Figure 5 . 10 : Heat of mixing for methanol + cyclohexane at T=416.29 K and 1900 kPa . 96 Figure 5 . 11 : Calculated and experimental activity coefficients for methano l + cyclohexane at 101.3 kPa . Data and fits correspond to the experimental VLE data in Figure 5 . 9 Figure 5 . 12 : Percent errors of K - ratios for methanol + cyclohexane system . Figures (a) and (c) are K - ratios in the VLE region and (b) and (d) are K - ratios in the LLE region (a) (b) (c) (d) 97 The methanol + n - pentane system is described in Figures 5 . 13 - 5 . 15 . The parameters were fitted to data at two temperatures, 372.7 and 397.7 K, and the results at both conditions were very similar. The higher temperature fits are shown here. It is clear that both CPA and PC - SAFT have difficulty describing the dew line for methanol + n - pentane. In contrast, WAG - Nagata and WAG - SH capture both phase lines very well with only two fitted parameters each. Figure 5 . 13 : Phase equilibria diagram for methanol + n - pentan e at T=397.7 K Figure 5 . 14 : Calculated and experimental [120] activity coefficients for methanol + n - pentane at T=397.7 K . Data and fits correspond to the e xperimental VLE data in Figure 5 . 13 . 98 Figure 5 . 15 : Percent errors of K - ratios for methanol + n - pentane system . Figure (a) is the K - ratio of methanol and (b) is the K - ratio of n - pentane in the VLE region at T=397.7 K Ethanol + cyclohexane fits are provided in Figures 5 . 16 - 5 . 19 . This system has no liquid - liquid immiscibility. Thre e models perform particularly poorly for this system: NRTL(2), NRTL(4) and CPA. The first two models falsely predict an LLE and the last miscalculates the azeotropic composition. WAG - SH and the four parameter WAG models describe the VLE and heat of mixing data much better for this system. (a) 99 Figure 5 . 16 : Phase equilibria diagram for ethanol + cyclohexane at 101.3 kPa . Figure 5 . 17 : Heat of mixing for ethanol + cyclohexane at T=298.14 K and 101.3 kPa . 100 Figure 5 . 18 : Calculated and experimental activity coefficients for ethanol + cyclohexane at 101.3 kPa . Data and fits correspond to the experimental VLE data in Figure 5 . 16 . Figure 5 . 19 : Percent errors of K - ratios for ethanol + cyclohexane system . Figures (a) is the K - ratio of ethanol and (b) is the K - ratio of cyclohexane in the VLE region at 101.3 kPa 101 Figures 5 . 20 - 5 . 22 show the model fits for ethanol + n - decane. The parameters were fitted to the data in the figures and two more P - xy VLE data sets in Narasigadu et al. [129] . The P - xy fits at the other temperatures are similar and, therefore, omitted for brevity. Linear segments in the bubble line or the heat of mixing indicate prediction of VLLE. Figure 5 . 20 : Phase equilibria diagram for ethanol + n - decane at T=338.17 K Figure 5 . 21 : Phase equilibria diagram for ethanol + n - decane at 101.3 kPa 102 Figure 5 . 22 : Heat of mixing for ethanol + n - decane at T=298.15 K and 101.3 kPa . The binary ethanol + n - decane system is especially challenging for all models tested. NRTL (4) incorr ectly predicts VLLE for the isothermal fits between 328 K and 348 K. WAG - NRTL (4) begins to incorrectly predict VLLE at a temperature slightly above 328 K and WAG - Nagata (4) does so at a temperature between 328 K and 298 .15 K. WAG - SH and PC - SAFT do not inc orrectly predict VLLE at 328 K but they cannot model the steep rise in the bubble line for the P - xy diagrams. The heat of mixing data shown in Figure 5 . 22 does not indicate phase splitting at 298.15 K, which would appear in the form of a linear section. Therefore, the VLLE predicted by almost all models tested here at temperatures above 298.15 K is incorrect and illustrates a shortcoming of the models. The CPA EOS generates a wider P - xy envelope than the other mode ls due to the difference in the association parameters used for CPA. The value used in CPA is approximately 1/5 th that of the value used in the WAG models and therefore, the hydrogen bonding is modeled to be weaker. The regressed residual portion must then compensate for the association 103 term in order to match the experimental data and the resulting residual contribution is overestimated. This demonstrates the importance of correctly quantifying the association effects on non - ideality and the diffic ulties that arise when the balance between and is incorrect. - NRTL and WAG - Nagata. For ethanol + n - decane, the experimental data used in the parameter regression was between 445 and 328 K. However, when the phase behavior was extrapolated 60 K lower, both models incorrectly predict a lower critical solution temperature (LCST) near 260 K. Interestingly, although CPA, PC - SAFT and WAG - SH provided inferior fits of the experimental data, t hey did not predict an LCST when extrapolated to low temperatures. We attribute this advantage to the temperature dependence of the residual contribution in the EOS models. Thus, future work will focus on studying the temperature dependence of the residual contribution for the local composition models. 5.3.3. Binary systems with self - and cross - associating systems Two systems with cross - association are modeled in this work. Figures 5 . 23 - 5 . 26 and 5 . 27 - 5 . 30 show the fits for methanol + water and methanol + ethanol respectively. For the first system, all the models describe the data well with no diffic ulty. However, the methanol + ethanol binary is more complicated because it is so close to an ideal solution. Almost all models can capture the phase behavior and heat of mixing for the system well except for CPA and WAG - SH. CPA is especially poor at model ing the system and Section 5.4 explores this limitation further. 104 Figure 5 . 23 : Phase equilibria diagram for methanol + water at 101.3 kPa Figure 5 . 24 : Phase equilibria diagram for methanol + water at T=298.14 K (bottom) and T=308.14 K (top) 105 Figure 5 . 25 : Calculated and experimental [121,122] activity coefficients for methanol + water at T=298.14 K, T=308.14 K and P=101.3 kPa from top to bottom . Data and fits correspond to the experimental VLE da ta in Figure 5 . 23 and Figure 5 . 24 . 106 Figure 5 . 26 : Percent errors of K - ratios for methanol + water system . Fi gures (a - c) are K - ratios of methanol and (d - f) are K - ratios of water in the VLE regions at T=298.14 K, T=308.14 K and P=101.3 kPa from top to bottom. (a) (b) (c) (d) (e) 107 Figure 5 . 27 : Phase equilibria diagram for methanol + ethanol at T=298.14 K (bottom) and T=413.13 K (top) 108 Figure 5 . 28 : Phase equilibria diagram for methanol + ethanol at T=298.14 K Figure 5 . 29 : Calculated a nd experimental activity coefficients for methanol + ethanol at T=298.14 K (bottom) and T=413.13 K (top) . Data and fits correspond to the experimental VLE data in Figure 5 . 27 . 109 Figure 5 . 30 : Percent errors of K - ratios for methanol + ethanol system . Figures (a) and (c) are K - ratios of methanol and (b) and (d) are K - ratios of ethanol in the VLE regions at T=298.14 K (top) and T=413.13 K (bottom) (a) (b) (c) (d) 110 5.3.4. Ternary sy stems The parameters that were obtained by regression to binary data are used to predict the behavior of the components in a ternary mixture. Mathias and Kister [134] have shown that including ternary LLE in the regression of binary parameters can lead to improved simultaneous representation. However, to explore the behavio r for predictive capabilities, the simultaneous regressions were not performed in this work. The phase diagram for methanol + n - heptane + cyclohexane is given in Figure 5 . 31 and the K - ratios for each of the compone nts are given in Figure 5 . 31 . The predicted LLE compositions are most accurate for the four parameter WAG - NRTL and NRTL models. WAG - SH and CPA also provide a good prediction. However, PC - SAFT and the other two - para meter WAG variations struggle to capture the methanol and cyclohexane rich phase. Figure 5 . 31 : Phase equilibria diagram for methanol + n - heptane + cyclohexane at 101.3 kPa and 298.15 K . Parameters used w ere fitted only to binary data. 111 Figure 5 . 32 : K - ratios for methanol (a) + n - heptane (b) + cyclohexane (c) system at 101.3 kPa and 298.15 K . 112 Similarly, the ternary diagram for methanol + ethanol + cycloh exane and K - ratios for the components is given in Figures 5 . 33 - 5 . 34 . The performance of the models is much different for tween ethanol and cyclohexane creates the wrong envelope in the figure. While this is the worst description of data, the other models are also weak at predicting the behavior of this system. However, WAG - SH provides the best prediction, out - performing CPA and PC - SAFT. Figure 5 . 33 : Phase equilibria diagram for methanol + ethanol + cyclohexane at 101.3 kPa and 298 K . Parameters used were fitted only to binary data. 113 Figure 5 . 34 : K - ratios for methanol (a) + ethanol (b) + cyclohexane (c) system at 101.3 kPa and 298 K . 114 5.3.5. General findings and recommendations Some general conclusions can be drawn from the results shown in this work. The first is that 4 - paramet er variations of WAG - NRTL and WAG - Nagata models provide excellent fits of experimental data for binary systems. The total SSQ of these two models is 47% and 15% smaller than that of 4 - parameter NRTL, which is the model with the next smallest SSQ, and are t herefore recommended as four parameter models. For systems that are difficult to model, the performance of WAG - NRTL model could be further improved by fitting the parameter. However, because the goal of this work is to compare WAG - Nagata and WAG - NRTL with the same number of adjustable parameters, this was not examined for the systems studied here, and is therefore left for future studies. Amongst two - parame ter models, the WAG - SH model developed here is a clear frontrunner. WAG - SH has been shown to correlate the behavior of the systems in this work better than the equations of state. It is particularly interesting to compare it to CPA. The functional forms of the energetic contributions in both WAG - SH and CPA are very similar but the SSQ for WAG - SH is only about 70% of that for CPA. The differences are that WAG - SH uses a 2/3 rd power on the umes are calculated Waals equation [8,67] . The parameter of is small for WAG - SH, CPA and PC - SAFT across all systems. However, the parameter of WAG - SH is consistently larger than that of CPA and PC - SAFT, in some cases by three orders of magnitude. This is because, the parameters in the EOS are corrections to the attractive parameter which have an imbedded temperature dependence in their functional form. In contrast, the fitted parameters in WAG - SH, and indeed all the WAG models, provide all of the 115 temperature dependence of the residual t erm. Therefore, the difference in the magnitude of the parameters between the WAG and EOS methods is expected. The small values of the fitted parameters in the equations of state indicate that the functional form of the equations adequately captures the te mperature dependence of the experimental data. However, the fits generated by them are still weaker than those by the 4 - parameter models for most systems. This is a consequence of the poor compositional dependence of CPA and PC - SAFT compared to WAG and NRT L. As indicated by Flemr [135] ) given by. 5 . 38 is the total number of molecules and is the number of molecules around . Lin et al. [136] showed that the non - com Boltzmann weighing factors for local compositions. Instead, they introduced a flexible weight, characteristic of each species, to enforce pair conservation. The resulting function, while more fundamentally correct, requires iterative solution. Moreover, similar to EOSs, the equation is less flexible with respect to composition than NRTL or N agata. Despite this, an interesting future with the combinatorial and association contributions used in this work. The SSQ analysis reported here is only s trictly applicable to the current set of systems and will therefore not necessarily represent all systems. However, the relative performances of the models are expected to yield similar trends. Further, some conclusions for this work may be dependent on th e association parameter. An observation in this work is that experimental alcohol + inert infinite dilution activity coefficients increase rapidly at infinite dilution. When the phase behavior is fitted using a model without chemical association, false LLE can be predicted because the residual 116 models do not include this rapid rise at infinite dilution (except for Nagata). However, adding chemical association has a strong effect at infinite dilution, and suppresses the tendency for false LLE because the asso ciation term creates very large infinite dilution activity coefficients before LLE begins to occur. Larger values of association constants permit fitting of infinite dilution activity coefficients while diminishing the tendency for false LLE and warrants f urther study. Future work will also investigate the relative magnitudes of the contributions to the activity coefficient, . For the WAG models and the systems tested here involving alcohol + hydrocarbons, the activity coefficient is dominated by the association contribution. The ratios of to the total at infinite dilution compositions are given in Table 5 . 8 . Table 5 . 8 : Association contribution to the total activity coefficient at infinite dilution System Component at infinite dilution WAG - N RTL(4) WAG - Nagata(4) WAG - SH WAG - NRTL(2) WAG - Nagata(2) Methanol + n - Heptane T=335 K Methanol 1.255 1.234 1.177 1.082 1.145 n - Heptane 1.368 1.274 1.069 1.070 0.903 Methanol + Cyclohexane T=335 K Methanol 0.988 0.958 0.954 - 0.924 Cyclohexane 0.880 0.7 93 0.725 - 0.712 Methanol + Water T=350 K Methanol 2.844 3.845 3.905 2.848 4.129 Water 1.226 1.793 1.881 2.055 1.794 Methanol + n - Pentane T=370 K Methanol - - 1.061 1.073 1.089 n - Pentane - - 0.897 1.060 0.955 Methanol + Ethanol T=340 K Methanol 1.0 56 1.048 0.935 1.025 1.028 Ethanol 1.064 1.060 0.902 1.029 1.030 Ethanol + Cyclohexane T=345 K Ethanol 0.877 0.861 0.851 0.773 0.840 Cyclohexane 0.680 0.684 0.711 0.725 0.641 Ethanol + n - Decane T=380 K Ethanol 1.429 1.448 1.245 1.211 1.212 n - Decan e 1.680 1.283 1.392 1.434 1.281 117 In general, the ratios are very similar across all WAG models for every system. However, there is significant variation in the calculated values for the methanol + water system. Since the calculated values of and are constant for a given system and condition, t his is an indication that the calculated energies of the dispersion in teractions, described by , are very different depending on the calculation method used. Compared to Nagata and SH, NRTL calculates larger deviations from ideality at infinite dilution due to dispersion interactions. Thus, future work should incl ude consideration of infinite dilution values. Association parameter value The parameters used to model ethanol with CPA were different than those used to calculate the ethanol association parameters in WAG. In the former case, the parameters, which will be referred to here as PM1, were calculated by fitting to vapor pressure and saturated liquid density data [1] . This procedure is typical of equations of state such as PC - SAFT and CPA. However, o ther studies have employed different methods for the same purpose. For example, Ren on and Prausnitz [68] and Campbell [16] fit the association parameters directly to binary phase equilibrium data. O thers [70,96] follow an approach developed by Brandani et al. [73] that calculates the association parameter by quantifying the difference between the vapor pressure of an alcohol and that of an ether homomorph with a similar molecular size. In a recent study, Kontogeorgis et al . [2] refitted CPA parameters to data that included the fraction of molecules that are non - bonded at equilibrium (monomers). This quantity was calculated from spectroscopy and the resulting parameter set, PM2, was found to be significantly different from those calculated previously. The values fitted to spectroscopy for ethanol give more associ ation and thus a larger association contribution to infinite dilution. Unfortunately, the parameters also do not simultaneously represent the vapor pressure. 118 We illustrate this difference for ethanol + cyclohexane and ethanol + n - decane in Figures 5 . 35 - 5 . 36 . Figure 5 . 35 : Comparison of phase equilibria fits with two sets of CPA parameters for ethanol + cyclohexane at 101.3 k Pa . Solid lines are fits calculated using PM1 parameters [1] and dashed lines are calculated with PM2 parameters [2] Figure 5 . 36 : Comparison of phase equilibria fits with two sets of CPA parameters for ethanol + n - decane at T=338.17 K . Solid lines are fits calculated using PM1 parameters [1] and dashed lines are calculated with PM2 parameters [2] 119 The solid and dashed lines in Figures 5 . 35 - 5 . 36 are fits of the CPA model made using PM1 and PM2 parameters respectively. The end points for pure ethanol are completely missed by the PM2 parameters. Thus, PM1 parameters were implemented for CPA fits in this work. Exploratory work show ed that most of the WAG models resulted in false LLE when regressed with PM1 for ethanol + cyclohexane. Also, because of the direct relationship between and monomer fraction, the PM2 parameters, which are fitted to the monomer fraction, are expected to be more fundamentally consistent with reality. Therefore, we elected to use PM2 parameters for ethanol to calculate the association parameter for WA G models. The parameters for the self - association of primary alcohols as calculated by the various approaches described are shown in Figure 5 . 37 . Figure 5 . 37 : for the self - association of the primary alcohols at 298.15 K calculated by fitting to different properties . Parameters were calculated by fitting to pure component properties (CPA [1,64] , PC - SAFT [115] ), binary VLE data (Renon [68] , Campbell [16] ) or vapor pressures of ether homomorphs (Bala [96] , using methodology from Nagata [70] ). The marker × spectroscopic and pure component data 120 For the methods of Renon, Campbell and Bala, association was modeled using chemical theory and the calculated parameter was the equilibrium constant, . In these cases, we use the finding of our earlier work [96] which showed that, with reasonable assumptions, , to compare the values of the . It is immediately clear from Figure 5 . 37 that the value of can vary by a factor of over 6 depending on the calculation method used. This is due, in part, to the fact that all these methods, with the exception of the one that yields the PM2 CPA parameter set, lack molecular - level experimental insight. Instead, they depend on macroscopic manifestations of hydrogen bonding. To address this situation, more fundamental techniques such as spectroscopy must be utilized in parameter calculations a nd we explore these methods in the following chapters. It is worth noting that calculated for the self - association of ethanol from the PM2 parameters, which include spectroscopic information, falls more in line with the values calculated from other methods and are close to that calculated by PC - SAFT. The work of Asprion et al . [7,137,138] uses s pectroscopy, but is not directly useful in the Wertheim framework. Regardless of the uncertainty in the exact value of the association constant, this work demonstrates that inclusion of association is important to improve the performance of traditional act ivity coefficient models. Conclusions In this work, a Wertheim activity coefficient term is combined with existing thermodynamic models including NRTL, Scatchard - resulting models are used to fit phase equilibria and heat of mixing data for methanol and ethanol - containing systems and compared to CPA and PC - SAFT. The four parameter WAG models were found to perform the best, with sum of squares errors that were 3 to 22% of those of 2 - parameter models . Among the two parameter models, WAG - SH provided comparable or, in some cases, 121 superior results to CPA and PC - fit simultaneous VLE and LLE data even though there was uncertainty in the associa tion constant value. While CPA parameters from literature were used to calculate the association parameter values in this work, the resulting fits highlight the importance of incorporating spectroscopic findings in determining accurate values of the associ ation parameter . 122 APPENDIX 123 Appendix H : s Table H. 1 : Regressed parameters values . The parameter was not regressed but is tabulated for conv enience Methanol ( i ) + n - heptane ( j ) WAG - NRTL(2prms) - - 530 - 318 0.2 WAG - NRTL(4prms) - 2.71 0.233 1400 - 426 0.2 WAG - SH - 0.0505 24.7 - CPA - 0.0298 0.0336 - WAG - Nagata(2prms) - - - 140 - 17 .5 - WAG - Nagata(4prms) 4.85 - 0.691 - 1532 177 - NRTL(2prms) 0 0 541 387 0.2 NRTL(4prms) 0.781 - 3.96 335 1515 0.2 PC - SAFT 0.0830 - 0.0455 - Methanol ( i ) + n - pentane ( j ) WAG - NRTL(2prms) - - 374 - 225 0.3 WAG - SH - 0.128 57.7 - CPA - 0.0186 0.0557 - WAG - Nagata(2prms) - - - 79.4 - 22.2 - NRTL(2prms) - - 552 492 0.3 PC - SAFT 0.0534 0.00770 - Methanol ( i ) + Cyclohexane ( j ) WAG - NRTL(4prms) - 0.645 - 0.643 732 - 101 0.2 WAG - SH - 0.0965 46.8 - CPA** 0.0320 0.00782 - WAG - Nagata(2prms) - - - 14.1 - 80.1 - WAG - Nagata(4prms) 1.62 - 0.196 - 502 - 23.7 - NRTL(2prms) - - 354 459 0.2 NRTL(4prms) 0.0803 - 4.780 347 1906 0.2 PC - SAFT 0.0726 - 0.0220 - Ethanol ( i ) + n - decane ( j ) WAG - NRTL(2prms) - - 41.2 - 20.4 0.2 WAG - NRTL(4prms) - 6.61 5.56 2170 - 1746 0.2 WAG - SH - 0.0154 7.10 - CPA - 0.122 0.122 - WAG - Nagata(2prms) - - - 5.63 - 18.1 - WAG - Nagata(4prms) 5.635 - 1.732 - 1721 472 - NRTL(2prms) - - 138 749 0.2 NRTL(4prms) 1.35 - 0.554 - 63.0 690 0.2 PC - SAFT 0.0634 - 0.0363 - 124 Table H. 1 Ethanol ( i ) + Cyclohexane ( j ) WAG - NRTL(2prms) - - 201 - 55.6 0.3 WAG - NRTL(4prms) 0.466 - 0.529 292 - 52.5 0.3 WAG - SH - 0.0189 17.4 - CPA - 0.0560 0.0947 - WAG - Nagat a(2prms) - - - 152 - 19.6 - WAG - Nagata(4prms) 0.3605 - 0.103 - 256 16.8 - NRTL(2prms) - - 187 794 0.3 NRTL(4prms) - 0.00483 0.912 238 415 0.3 PC - SAFT 0.0646 - 0.0225 - n - Heptane ( i ) + cyclohexane ( j ) WAG - NRTL(2prms)* 0 0 - 95.1 123 0.3 WAG - NRTL(4prms) 0.0367 - 0.494 - 142 344 0.3 WAG - SH - 0.0202 8.41 - CPA 0.0123 - 0.0119 - WAG - Nagata(2prms)* 0 0 - 193 184 - WAG - Nagata(4prms) - 0.159 0.625 65.4 - 259 - NRTL(2prms)* 0 0 - 123 149 0.3 NRTL(4prms) - 0.0382 - 0.435 - 131 334 0.3 PC - SAFT - 0.00573 0.00787 - Methanol ( i ) + Water ( j ) WAG - NRTL(2prms) - - - 473 322 0.2 WAG - NRTL(4prms) 12.7 - 4.74 - 2881 720 0.2 WAG - SH - 0.147 - 26.6 - CPA - 0.0881 0.0163 - WAG - Nagata(2prms) - - 130 59.6 - WAG - Nagata(4prms) 1.72 - 2.76 - 489 1058 - NRTL(2prms) - - - 111 345 0.2 NRTL(4prms) 3.59 - 3.53 - 1224 1430 0.2 PC - SAFT 0.00387 - 0.0759 - Methanol ( i ) + Ethanol ( j ) WAG - NRTL(2prms) - - 67.8 - 63.2 0.3 WAG - NRTL(4prms) - 0.868 0.712 336 - 291 0.3 WAG - SH 0.0145 1.28 - CPA** 0.0155 - 0.0365 - WAG - Nagata(2prms) - - - 49 .4 32.2 - WAG - Nagata(4prms) 1.03 - 0.660 - 327 215 - NRTL(2prms) - - 62.8 - 53.8 0.3 NRTL(4prms) - 1.03 0.929 351 - 320 0.3 PC - SAFT 0.00926 - 0.00830 - 125 Table H. 2 : Hayden - parameters . The value of , which is dimensionless, for alcohol + alkane binaries is zero Comp j Methanol Ethanol Comp i Methanol 1.63 1.55 Ethanol 1.55 1.4 126 Chapter 6. Quantitative Analysis of Infrared Spectra of n - Butanol + Cyclohexane with Quantum Chemical Calculations Hyd rogen bonding has profound effects on the microscopic behavior of molecules. Infrared spectroscopy (IR) allows for the analysis of molecular vibrations by excitation with infrared radiation. However, quantitative analysis of the hydroxyl band in the IR spe ctrum, where hydrogen bonding is most prominently expressed, is non - trivial. Specifically, the broadness of the band and the range of variation of the O - H absorption coefficient complicate the analysis. In the present work, sequential MD and QM simulations are used to develop functions that relate the vibrational frequencies to the integrated absorption coefficients of hydroxyl bands. This relationship is then used to quantitatively calculate the mixture concentration of alcohol molecules from experimental IR spectra by integration across the entire hydroxyl band. Finally, universal maps to relate bond lengths, vibrations and NMR chemical shifts are proposed. 6.1. Introduction In infrared (IR) spectroscopy, molecules are excited with light in the mid IR wavelen gth region and the amount of light absorbed by the sample is detected. Absorption of light at different wavelengths induces vibrations of different bonds within the molecule and analysis of the resulting spectra can yield information about the types of fun ctional groups it possesses. IR spectroscopy has long been used in this way to identify chemical compounds through the interpretation of band patterns. Moreover, quantitative analyses of IR spectra can be used to gain insight into the concentration of func tional groups in a solution. For example, Williams et al. [139] explored the relationship between ab solute integrated intensities of the C - H stretching and bending bands of gas - phase alkanes. The authors compared density functional theory (DFT) calculation results to experimental IR spectra and found that the number of C - H bonds in the molecules studied is 127 linearly correlated with the integrated intensities of C - H stretching and bending bands. IR is also frequently used in chemometrics [140] . Several researchers have been able to correlate intensities of various classes of compounds with physical characteristics of the molecules such as the number o f methylene groups [141] , molecular size [142] and degree of branching [143] . The O - H stretching bands have been the focus of many studies that aim to understand the complex effects of hydrogen bonding [37,144 151] . These bands occur between ~3200 cm - 1 and 3700 cm - 1 and their combinations appear in the form of two overall peaks: a sharp higher frequency peak and a broad lower frequency peak. The formation of hydrogen bonds is known to decrease the O - H stretching band frequency and increase its integrated intensity [152] . Wu et al. [146] investigated supercritical and liquid methanol and found that an isothermal increase in density causes the integrate d spectroscopic area to increase and the hydroxyl vibrational band to shift to lower frequencies. Expectedly, isobaric heating has the opposite effect. Due to these widely varying measurement response factors compared to C - H bonds, vibrations associated wi th free and hydrogen - bonded hydroxyl groups are much more difficult to interpret in quantitative terms and have therefore largely been interpreted qualitatively. To overcome some of the challenges involved in IR analysis of hydroxyl peaks, computational to ols such as molecular dynamics (MD) [153] and quantum mechanical (QM) [139,144,154 159] simulations are being leveraged to elucidate the effects of hydrogen bonding on IR peak characteristics. MD and QM calculations balance computational expense with modeling rigor. The former method is less deman ding of computational time and resources but does not model electronic effects, which are the basis of hydrogen bonding and, in particular, result in the distributions of hydrogen bond types. Indeed, Kwac and Geva [17] show that, depending on the choice of empirical forcefield, the simulated 128 relative populations of species can vary by as much as 20%. For this reason, MD is often used in conjunction with QM [17,160 162] rather than as the primary tool for investigation. A recent development in the computational sphere is empirical mapping. This approach, developed by Skinner et a l. [163 167] frequencies to properties that can be calculated with reasonable accuracy from MD calculat ions. without having to use excessive computational resources. Mesele and Thompson [168] extended frequencies, dipole derivatives and position matrix elements to the electric field on the atoms. In this work, we aim to present a combin ed computational and experimental approach which leverages the power of simulations to address the challenges of interpreting the infrared behavior of O - H bonds. Furthermore, we apply this technique to quantitatively analyze the entire hydroxyl band and ca lculate the relative and absolute concentrations of hydroxyl groups in the various contexts (monomers and clusters) existing in solution. In all our discussions, references to the hydroxyl vibrational bands pertain to the vibrations of the covalent O - H bon d and not vibrations of the actual hydrogen bond (which appear at much lower frequencies [152] ). Additionally, the gnoring the clusters formed by hydrogen bonding the molecules are counted individually. The term chemistry literature. 6.1.1. Analysis of the O - H stretching band Qua ntitative interpretation of infrared spectra begins with the Beer - Lambert law, which relates the observed IR absorbance of a solute to its concentration in solution according to: 129 where and are the observed absorbance and concentration of the absorbing solute i respectively. The pathlength, , is the length of the sample that t he IR light beam passes through and is the molar absorption coefficient, als o known as the extinction coefficient, of species i . The latter parameter is a measure of the IR absorbance of a single bond vibration at a selected wavenumber and is often assumed to be independent of concentration for the solute. In common practice, solu tions of known concentrations are prepared and analyzed with IR spectroscopy. Then, a peak corresponding to a vibration in the solute molecules is chosen in an area where overlap with solvent bands is minimal. The absorbance peak height values are then plo tted against the experimental concentrations of the solutions. Finally, the value of is calculated as the gradient of this plot, which is ideally linear, and used in subsequent studies to analyze solutions of unknown concentration. However, the Beer - La mbert law is not universal; it fails under conditions where the relationship between and becomes non - linear. This occurs when the solute molecules interact significantly with one another or solvent molecules, such as in high concentration s olutions. In these cases, the value of the molar absorption coefficient, , can vary dramatically for the same bond depending on the molecule conformation and its microenvironment. This contributes to the broadening of the observed IR band and complicates quantification. Moreover, there is disagreement in the literature concerning the assignment of vibrational bands. In earlier studies [137,145,169] , vibrational bands were assigned to hydrogen bonded clusters and were distinguished based on the size of the cluster. Hall and Wood [170] were the first to propose that covalent O - H bond vibrations should be classified individually according to if and how t hey participate in hydrogen bonding. If the O - H is neither accepting nor donating a h ydrogen atom, it 6 . 1 130 Figure 6 . 1 . Using femto s econd IR pump probe studies of ethanol clusters, Woutersen et al. [171] were able to assign t he , , and bonds to experimentally deconvoluted peaks as far back as 1997. Figure 6 . 1 : Types of covalent O - H bonds This approach to IR band interpretation is more aligned with the funda mental basis of IR excitation as a bond perturbation rather than a molecular or cluster - wide phenomenon. Therefore, Hall and [147,148,150,172] and will also be used in this work. 6.1.2. Motivation While research in this area is extensive, we are unaware of any work demonstrating quantification of the entire hydroxyl IR band area for alcohols and relating the area to the apparent concentration. Such a relationship has profound implications for the s tudy and modeling of hydrogen bonding. Indeed, in the development of engineering models for the association of an alcohol in an inert solvent, the key parameter that defines the extent of hydrogen bonding is the fraction of hydroxyl protons that remains no nbonded at equilibrium. Once this value is determined, some assumptions Accepts: 0 Donates: 0 Accepts: 0 Accepts: 1 Ac cepts: 1 Donates: 1 Donates: 1 Donates: 0 Accepts: 2 Donates: 1 Accepts: 2 Donates: 0 131 can be made about the types of species present in solution and the energy of their formations. A simple model that has often been adopted [98,173,174] for alcohol + inert systems is the CLAM, or Continuous Linear Association Model. As the name suggests, CLAM accounts only for linear species in solution and assumes that each hydrogen bond has the same energy of formation regardless of the size of the cluster it forms. The motivation be hind the current study is to develop a procedure capable of accurately determining the fraction of free hydroxyl protons and conduct a thorough quantitative analysis of the O - H IR bands. 6.2. Computational m ethods To analyze the vibrational spectra derived fr om infrared spectroscopy, a combined molecular dynamics (MD) and quantum mechanics (QM) approach was used. In this section, these computational techniques are outlined. 6.2.1. Molecular d ynamics s imulations Molecular dynamics simulations were carried out using t he AMBER 14 package [175] . The AMBER94 force field was implemented with the AM1 - BCC charge method and no modification to the force field parameters. For each concentration, a cubic box of n - butanol and cyclohexane molecules was created at the ideal solution density using PACKMOL [176] and the energy of the system was minimized within 1500 steps. Next, the box was heated up with a 40 ps NVT run, using a time step of 2 fs, in two stages. The temperature was ramped up from 0 to 283.15 K during the first 9000 steps then maintained at that temperature for the remaining 11000 steps. Finally, a 10.4 ns NPT production run was conducted at 1 b ar and a time step size of 2 fs. The temperature and pressure were controlled using the Langevin thermostat (with a collision frequency of 2 fs) and Berendsen barostat respectively; both implemented with default parameters. During the NVT and NPT simulatio ns, periodic boundary conditions were enforced in the x,y and z coordinates. 132 The cutoff for non - bonded interactions was set at 8 Å and a continuum model was used for Lennard - Jones interactions beyond this range. Electrostatic interactions were calculated u sing the Particle Mesh Ewald (PME) summation method with default settings and parameters. Only bond lengths involving hydrogen were constrained using the SHAKE algorithm. This method was repeated for n - butanol + cyclohexane and ethanol + cyclohexane at the equimolar composition. The simulation details that varied between systems are given in Table 6 . 1 . Table 6 . 1 : MD simulation d etails . Box and run for n - butanol + cyclohexane systems studied in this work System n - Butanol + Cyclohexane Ethanol + Cyclohexane Alcohol mole fraction 0.1 0.5 0.5 Number of alcohol molecules 26 128 168 Number of cyclohexane molecules 230 128 168 Ini tial cubic box length (Å) 40 40 36 Temperature (K) 283.15 283.15 298.15 Production run length (ns) 10.4 10.4 13.0 T he production runs were longer than 10 ns in all cases, giving all the systems ample time to equilibrate. Furthermore, only the last 2.4 ps of the simulations were analyzed for hydrogen bonding information. The hydrogen bond criteria were defined as an O - O distance < 3.2 Å and an O - nd each bond was assigned a class based on the categories in Figure 6 . 1 . For the current study, only linear species are analyzed. Furthermore, the low incidence of structures including and bonds suggested tha t they play a minor role; they are thus neglected. To verify thermal equilibration, Figure 6 . 2 shows the number of times each n - butanol molecule appears as a , or bond in trimers in approximately 590 evenly spaced frames collected every 40 fs in the last 2.4 ps of the production run. Each molecule appears in a trimer in about 15% of 133 the frames. The x - axis is the unique identifying number of each n - butanol molecule in the box. The distribution is random for t he three bond types across all molecules in solution. Figure 6 . 2 : Bond distribution in trimers calculated with MD simulations for an equimolar n - butanol + cyclohexane mixture . The distribution illustrates that molecules are equally likely to sample all bond types. 6.2.2. Quantum m echanics s imulations Once the hydroxyl bonds were classified, each cluster was prepared to undergo QM calculations. Careful implementation of QM simulations is a balance between accurac y and computational expense. For this reason, performing extensive calculations over the entire box was not feasible and an alternative was necessary. In this work, all molecules with atoms that fell within 5 Å of the hydroxyl hydrogen atom in the hydrogen bonded cluster were retained in the simulation environment. All other molecules were excluded. This cutoff effectively ensures that the sampled environments were representative of the whole box while minimizing the computational effort 134 required for quantu m calculations. To this end, we also constrained all atoms and bonds except for the - CH 2 - O - H groups in each alcohol molecule of the hydrogen bonded cluster of interest. Gaussian 09 [177] was used to optimize the geometries of each cluster and perform the frequency calculations. The B3LYP [178] level of theory and 6 - 31G* basis set were chosen shown [154,159] to capture the effects of hydrogen bonding accurately for a reasonable computational cost. Finally, the IR vibrational information calculated underwent two tests. In the first, the final atom positions were checked to ensure that the covalent bond classifications made e arlier had not changed during QM optimization. Secondly, vibrational coupling between multiple - O - H bonds was estimated by calculating the displacement of each - O - H hydrogen atom in the cluster using: 6 . 2 where , and is the displacement of the hydrogen atom in the , , and directions respectively. Several studies have shown that vibrational coupling between different groups have an effect on the spectra [149,172,179,180] . However, in this work, we are interested in comparing the behaviors of individual bond types. Therefore, to reduce the effects of coupling, vibrations were excluded from analyses if they could not be confidently assigned to only one O - H bond. a vibration being at least 0.3 Å 2 6.3. Experimental m ethods Experimental IR data collection was performed by William G. Killian, a member of the Lira Lab at Michigan State University and will be included in h is Ph . D . dissertation [181] . Therefore, this section provides a basic outline of the experimental methodology for context. Cyclohexane, ethanol and n - butanol were purchased in anhydrous forms at purities of 99.5%, 99.5% and 99.8% respectively and further dried using 3 Å molecular sieves for at least 72 hours. Samples were 135 prepared volumetrically in a glove box and then transferred to scintillation vials to be removed from the inert environment for analysis. Infrare d absorbance spectra were collected on a JASCO FT/IR - 6600 spectrometer using a Specac demountable temperature - controlled liquid flow cell (GS20582) with CaF 2 windows and PTFE spacers. The system was allowed to stabilize for 10 minutes at each temperature. The pathlength for all measurements was 0.01092 cm. A background of the empty cell was taken at each temperature and subtracted from the sample data. For each reported spectrum, 128 scans were accumulated with a resolution of 2 cm - 1 . 6.4. Results and d iscussi ons 6.4.1. Processing and p reliminary a nalysis of e xperimental IR s pectra The raw IR spectra were processed as follows. The hydroxyl band region is determined to be 3049.9 to 3755.2 cm - 1 . First, the solvent bands were removed from the spectra by subtracting conc entration - weighted spectroscopic data of neat cyclohexane at the same temperature and pathlength of the sample. The concentration was calculated by assuming ideal volume of mixing using temperature - dependent volumes through . The pro cessed experimental IR spectra for n - butanol in cyclohexane in the region of the hydroxyl stretching band is given in Figure 6 . 3 for = 0.1 . As the temperature increases, the peak at ~3650 cm - 1 and sho ulder at ~3540 cm - 1 , and respectively, increase in absorbance while the broad peak centered at ~3320 cm - 1 , , diminishes. 136 Figure 6 . 3 : Experimental O - H IR band for a 0.1 mole fracti on of n - butanol in cyclohexane . Data was collected at five temperatures between 30 C and 70 C While band assignment is a topic of great interest to the scientific community, we defer such discussion to the following chapter. In this work, we instead co nduct a quantitative analysis of the entire hydroxyl band. To this end, we begin with an investigation of the physical significance of the absorbance band area. Figure 6 . 4 shows the relationship between the area u nder the entire O - H absorbance band and the apparent concentration for four concentrations at five temperatures. 137 Figure 6 . 4 : Total absorbance band area of the O - H band for n - butanol + cyclohexane data as a function of the mole fraction for 4 compositions . Spectra were collected at = 0.0469,0.0678, 0.0817 and 0.100 and the lines indicate the line of best fit given by the equations on the figure. A thermal correction for density is applied. Whe n temperature rises, it is expected that an increasing number of hydrogen bonds will break, forming smaller chains and monomers. Thus, the distribution of the different bond types will change. However, at a given composition, the total quantity of hydroxyl groups should remain constant across the temperature range studied here. To correct concentration for the effect of temperature on the density, the band areas plotted in Figure 6 . 4 are scaled with a thermal correc tion factor which is defined as: 6 . 3 where and are the mixture densities at the temperature of the experiment and a reference temperature (298.15 K) respectively. The m ixture density is calculated by assuming no excess volume. Therefore, where is the molar volume of component i . y = 205.41x - 49.18 y = 220.56x - 40.32 y = 235.52x - 31.51 y = 253.65x - 25.47 y = 269.34x - 19.34 0 50 100 150 200 250 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 Total OH Absorbance Band Area (cm - 1 ) Mole fraction of n - butanol 138 It is immediately clear from Figure 6 . 4 that the assumed linear ity of and that is traditionally assumed for a given absorbance in the Beer - Lambert law does not apply for the overall hydroxyl band. If it did, all of the points (at all temperatures and compositions) would be reasonably well represented w ith a single straight line that passes through zero on both axes. However, while data at each temperature is quite linear, the intercept varies significantly between temperatures, indicating that the changing distribution of bond types is not being capture d. This discrepancy is expected as hydrogen bonding is known to conflict with the assumptions that underlie the Beer - Lambert law. Figure 6 . 4 is evidence that the molar absorption coefficient, , must vary at different vibrational frequencies. Next, the results of the QM/MM analysis are shown to provide key insights into this variability, resulting in relationships between absorption intensity, wavenumber, and covalent bond length. 6.4.2. Results from MD+QM In this section, we investigate the calculated IR characteristics of the hydroxyl stretch. Because of the computational expenses associated with each QM calculation, the scope is limited to monomers, dimers, trimers, and tetramers and Table 6 . 2 lists the numbers of each bond and species type analyzed. Table 6 . 2 : Number of each species and bond analyzed with QM calculations . All species larger than a monomer hav e one and one bond. Trimers and tetramers also possess one and two bonds respectively. Cyclic clusters are not considered. 0.1 0.5 Total Species Monomer 653 366 1019 Dimer 153 230 383 Trimer 53 52 105 Tetramer 46 23 69 Bonds 653 366 1019 252 305 557 145 98 243 252 305 557 139 the integrated absorption coefficient, , in units of km/mol, which is defined by the equation: 6 . 4 In this function, the integral spans the entire absorption band. Like the molar absorption coefficient, the integrated absorption coefficient is directly related to the derivative of the dipole mome nt and is a measure of the intensity of a transition. Literature on the relationship between the integrated absorption coefficient, , and the molar absorption coefficient, , is somewhat scarce. However, Spanget - Larson [182] provides a simple interpretation of the process which involves assuming a function for . When a Loren tzian function is assumed for , then is related to the peak height selected for via a constant. F or the purposes of this work, the results will be scaled by a constant and thus avoiding the need for numerical conversion or specificat ion of the Lorentzian parameter values. We use IR absorbance, with the understanding that it is proportional to the peak height . Following standard practice, the vibrational frequencies resulting from Gaussian are scaled by a factor of 0.96 to correct for the limitations of the B3LYP/6 - 31G* model and better match experimental IR data [183] . To visualize general trends, the average vibrational wavenumbers and integrated absorption coefficients for different bonds and species are given in Table 6 . 3 . The first observation is that with an increase in the apparent concentration of alcohol in the simulated box, there is a slight red - shift in almost all the vibrational frequencies, most prominently those of the and bonds. In general, there is also a red - shifting tendency for bond vibrations as the size of the chain they are on grows. Indeed, the difference is most prominent between dimers and larger oligomers with dimer and bonds having vibrational frequencies that are, on average, ~9 and ~50 cm - 1 greater 140 than those of trime rs and tetramers respectively. Table 6 . 3 follows a trend that increases in the order though the differences between species are more complex. Table 6 . 3 : Average calculated vibrational frequencies and integrated absorption coefficients Species Bond (cm - 1 ) (km/mol) 0.1 0.5 0.1 0.5 Monomer 3602 3587 31.90 53.97 Dim er 3601 3583 52.33 72.32 3479 3471 376.5 411.0 Trimer 3589 3572 68.31 96.21 3418 3414 464.1 525.2 3410 3412 583.9 537.0 Tetramer 3586 3574 68.39 87.65 3436 3426 426.8 428.1 3387 3363 628.8 675.9 Figure 6 . 5 shows the vibrational frequency and integrated absorption coefficient associated with each hydroxyl vibration at both apparent concentrations, and . When the two concentrations are considered separately, we observe little difference in the trends, indicating that apparent concentration has little effect on the - relationship. Therefore, they are plotted together in Figure 6 . 5 and the different m arkers denote the four bond types studied here: (diamonds), (triangles), (squares), and (circles). 141 Figure 6 . 5 : Hydroxyl stretching frequencies and integrated absorption coefficients for two con centrations of n - butanol + cyclohexane mixtures calculated from QM simulations . Diamond, triangle, square and circle markers denote , , , and bond vibrations respectively. The solid line is a line of best fit determined visually. It is evident t hat the two bonds in which the hydrogen atom is free, and , overlap completely and are responsible for the sharp higher frequency peak in the O - H stretching band. Together, these two bonds constitute all the free hydroxyl protons in solution. This is convenient since the key parameter required for the thermodynamic modeling of alcohol + inert systems is the fraction of the free O - H protons . Moreover, there is significant disagreement in the literature concerning the - overlap in IR spectra. Several authors have assumed that bond vibrations do not contribute significantly to the sharp free O - H proton peak, instead allocating it entirely/predominantly to the vibration [184 186] . In these cases, it is assumed that most clusters in solution are in cyclic form (resulting in few bonds), that the peak occurs at different frequency altogether, or that the intensity of the absorbance is significantly less than that of the . However, the results of o ur QM calculations are consistent with more recent work in this area [2,150] . While there is a red - 0 200 400 600 800 1000 1200 1400 3200 3300 3400 3500 3600 3700 Integrated Absorption Coefficient (km/mol) Wavenumber (cm - 1 ) 142 shift in the vibration due to charge redistribution when an oxygen atom accepts another proton, this shift is very slight. This indicates that the electronic structure is only weakly affected when an oxygen in a hydrox yl hydrogen bonds. The next observation is that the bonds appear at a lower frequency and have greater integrated absorption coefficients, (approximately 3 - 10 times that of the and bonds). The bonds follow the same pattern with respect to bo nds with integrated absorption coefficients that are 3 - 20 times that of the and bonds). The frequency trend is easily bonding, the covalent bond is weakened. As a result, the potential energy well in which the hydrogen in moving is effectively broadened, causing the vibrational frequency of the hydroxyl to red shift. The increase in intensity is caused by an increase of the dipole moment due to charge distribut ion that occurs through the hydrogen bond as the proton vibrates between the two oxygen atoms. The most valuable finding of this work and, in fact, the most obvious is that the relationship between the integrated absorption coefficient of an O - H bond and its vibrational wavenumber follows a curve that is independent of the bond category. The relationship can be described well by an inverted Gaussian function of the form: 6 . 5 where the parameters , , and are visually determined. This fit is plotted as a black line in Figure 6 . 5 . It is worth noting t hat a polynomial function was insufficient in capturing the curve characteristics well. Previous studies in this area have recognized patterns in vibrational characteristics. For example, as early as 1956, Huggins and Pimentel [187] , using notation X - Y to indicate that a variety of electronegative groups, found similarly interesting patterns for a 143 wide range of hydrogen - bonding systems. They show that a monotonic relationship between the vibrational X - H frequency and , which is defined as the difference between the intensity in the solvent minus the intensity in chloroform, a weak hydrogen bond donor. More recently, Mesele and Thompson [168] conducted DFT calculations on neat alcohols and showed th at the empirical maps that relate transition frequencies, position matrix elements, dipole derivatives and the electric field are surprisingly linear. Moreover, these relationships were identical for all 4 primary alcohols tested and were predicted to hold to all other alcohols. Having independently uncovered this pattern, we further explore the applicability of Figure 6 . 5 and Eq. ( 6 . 5 ) by repeating the described procedure for an equimolar mixture of ethanol and cyclohexane at 298.15 K. The resulting plot of vibrational frequencies and integrated extinction coefficient overlaps completely with data in Figure 6 . 5 as shown in Figure 6 . 6 for all bond types. Further, we also include the results of Murdoch et al [144] in which the vibrational behavior of ethanol clusters from monomers to hexamers was calculated wit h DFT methods. The clusters were optimized in vacuum and the effects of conformation (anti and gauche) were also investigated. In Figure 6 . 6 , the linear species are included as filled squares. For the cyclic specie s, we noted that, while there was a wide scatter in the intensities of the bond vibrations, each high intensity vibration was coupled with a low intensity vibration. From an experimental perspective, an IR experiment will collect the average intensity. Thu s, we averaged the values of frequency and integrated extinction coefficient for all the molecules in each cluster and included them in Figure 6 . 6 as filled triangles. The fit of these calculations with our own is striking. Though untested thus far, Figure 6 . 6 is a promising indication that the inverted Gaussian fit may be representative for other alcohols. 144 Figure 6 . 6 : Hydrox yl stretching frequencies and integrated absorption coefficients for other systems calculated from QM simulations . Unfilled diamond, triangle, square and circle markers denote , , , and bond vibrations of ethanol + cyclohexane respectively (this work). Filled small circles are butanol + cyclohexane data (this work). Filled squares and triangles denote linear and averaged cyclic ethanol oligomers respectively (Murdoch et al. [144] ). The solid line is a line of best fit determined visually. 6.4.3. Scaling and f urther a nalysis of e xperimental IR s pectra To quantify the number of bonds in the hydroxyl band in Figure 6 . 3 , we revisit the Beer - Lambert law and use the integrated extinction coefficient, in place of . At the wavenumber of every data point, a value for the is calculated and used to scale the absorbance. In addition to this, the wav enumber quantities found using QM calculations must be applied carefully when interpreting experimental data. It should be recognized that relative relationships between ab initio calculated IR features are more reliable than absolute values. Therefore. we adjust the inverted Gaussian function developed in Eq. ( 6 . 5 ) such that the wavenumber of the minimum, which is decided by parameter , matches the free - end peak center in the experimental spectra. Therefore, we set to 3645 cm - 1 and is calculated according to: 0 200 400 600 800 1000 1200 1400 3200 3300 3400 3500 3600 3700 Integrated Absorption Coefficient (km/mol) Wavenumber (cm - 1 ) 145 6 . 6 where represents the scaled absorbance at wavenumber and represents a concentration response by the system to IR irradiation. Plotting the resulting scaled spectra yields Figure 6 . 7 which is visually very different from the original spectra in Figure 6 . 3 . Peaks and now dominate in intensity while is diminished considerably. Figure 6 . 7 : Scaled absorbance spectra for n - butanol + cyclohexane at = 0.1 . Spectra are scaled w ith the integrated absorption coefficient function Calculating the area under the new scaled spectra and comparing the result s with the experimental mole fraction provides Figure 6 . 8 , which is analogous to Figure 6 . 4 . 146 Figure 6 . 8 : Scaled absorbanc e band area of the O - H band for n - butanol + cyclohexane data as a function of the mole fraction of 4 experimental compositions . Spectra were collected at = 0.0469,0.0678, 0.0817 and 0.100 and the line indicate the line of best fit given by the equations on the figure. A thermal correction for density is applied. The relationship between the band area and the apparent mole fraction is corrected significantly. First, at a given composition (and correcting for temperature effects on density due to thermal expansion), the integrated concentration response of the system is remarkably close regardless of temperature for all the compositions studied. T his is especially clear when the sample standard deviation between temperatures is compared as in Table 6 . 4 . The standard deviation is calculated using the equation: 6 . 7 where s is 5, the number of temperature data points at each composition. y = 368.19x R² = 0.9667 10 15 20 25 30 35 40 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 Scaled IR Band Area Mole fraction of n - butanol 147 Table 6 . 4 : Areas and standard deviations for the absorbance and scaled absorbance bands T( C) Absorbance Band, Scaled Absorbance Band, Area Area Std Dev Area Area Std Dev 0.0469 30 94.94 22.35 17.73 1.0013 40 82.05 17.74 50 68.35 17.54 60 5 3.85 16.98 70 38.40 15.37 0.0678 30 154.1 27.91 26.78 0.4453 40 138.1 27.08 50 120.2 26.82 60 102.0 26.52 70 83.99 25.91 0.0817 30 187.1 32.48 30.69 0.5723 40 168.9 31.11 50 149.2 30.87 60 126.7 30.16 70 105.6 29.69 0.1 30 229.9 34.68 35.04 0.3828 40 209.2 35.37 50 186.3 35.27 60 165.1 35.97 70 142.3 35.80 Figure 6 . 8 demonstrates that, with the scaled spectra, the Beer - Lambert law can be applied successfully to the entire O - H band permitting quantification of the hydroxyl band to determine the alcohol concentration in solution. The functional form of the inverted Gaussian function developed in Eq. ( 6 . 5 ) results adequate scaling. To our knowledge, no other work has performed a comparable quantitative analysis of the O - H band due to its known complexity. The work presented here provides a powerful tool that maps the band characteristics to the concentration in solution, effectively pla cing the absorptions of all O - H sites, regardless of their context, on a quantitatively equal footing. 148 6.4.4. Universality of s pectroscopic c haracteristics An interesting and important universality emerges from the QM calculations of various properties demonstra ted here for an equimolar mixture of ethanol + cyclohexane at T = 298.15 K. The same computational technique described in Section 6.2 were carried out on the system. Figure 6 . 9 shows ho w the length of the optimized O - H bond length correlates with its vibrational frequency. Figure 6 . 9 : Relationship between the covalent O - H bond length and vibrational frequency calculated for equimolar etha nol + cyclohexane mixture . Diamond, triangle, square and circle markers denote , , , and bond vibrations of ethanol + cyclohexane respectively The relationship is described extremely well by a linear function of the form: where is the length of the hydroxyl bond in Angstroms and is the wavenumber in cm - 1 . Furthermore, the IR vibration also shows a distinct relationship with the NMR chemical shift of the hydroxyl calculated with Gaussian09 as shown in Figure 6 . 10 . 0.96 0.965 0.97 0.975 0.98 0.985 0.99 0.995 3200 3300 3400 3500 3600 3700 Covalent OH Bond Length (Å) Wavenumber (cm - 1 ) 149 Figure 6 . 10 : Comparison of calculated NMR and IR characteristics for equimolar ethanol + cyclohexane mixture . Diamond, triangle, square and circle markers denote , , , and bond vibrations of ethanol + cyclohexane re spectively While some relationship between these different molecular characteristics is expected, the direct correlation is striking. Hydrogen bonding is often considered to be a complex phenomenon with random effects on mixture properties. These figures show that, with careful implementation and analysis, computational simulations have the power to create models that relate the bond length, of these elements from the others and can have far - reaching impacts on the analysis of the molecular behavior of alcohols. The visual fit of Eq. 6 . 5 may require refinement as more data become available at lower concentrations, or with other alcohols. T ypically, ab initio calculations at the level of theory implemented here would be expected to require additional scaling beyond the frequency shift that is implemented. However, with the currently available data, such scaling is difficult to justify -1 0 1 2 3 4 5 6 7 3200 3300 3400 3500 3600 3700 NMR chemical shift (ppm) Wavenumber (cm - 1 ) 150 becaus e the current constants provide scaled areas with excellent linear correlation to experimental data. In future work, we aim to leverage these maps to gain quantitative insight into the speciation that occurs due to hydrogen bonding. It will be important to study the behavior of cyclic species and include them in model development. Our preliminary work with Murdoch [144] species in Figure 6 . 6 indicates that bonds in cyclic clusters behave very similarly to bonds. For this reason, we expect that they will also be well - represented with the QM correlation developed here. However, the behavior of and bonds is less easy to anticipate without further calculations. Therefore, future directions will explore the infrared responses of these bonds. Finally, we aim to calculate the values of association parameters for a thermodynamic model developed in an earlier work [96] from spectroscopic findings. 6.5. Conclusions I n this work, the hydroxyl vibrational band is analyzed quantitatively for alcohol + inert mixtures using insights from QM simulations. First, MD simulations were carried out to generate sample environments around O - H bonds. Then, selected clusters were fur ther analyzed with QM calculations to gain an understanding of their spectroscopic characteristics. It was found that the integrated absorption coefficient and vibrational frequency of a group can be related by an inverted Gaussian function. Through a vari ation of the Beer - Lambert law, this function can be used to calculate the concentration of hydroxyl groups in solution from IR spectra. Finally, universal - H bond length, its vibrational frequency, integrated absorption coefficient and NMR chemical shift were developed. 151 Chapter 7. Integration of Quantum Calculations and Spectroscopy for Wertheim Alcohol Association In this work, we provide preliminary calculations and results for the derivation of Wertheim association parameters from spectrosco pic findings. IR and NMR spectra for mixtures of n - butanol + cyclohexane mixtures were collected and analyzed using a combination of thermodynamic modeling and insight from MD+QM calculations. While some results are discussed, the challenges encountered he rein lie beyond the scope of this dissertation. As such, recommendations are given for future directions that will help advance the current work and address its limitations. 7.1. Introduction 7.1.1. Hydrogen bonding and infrared spectroscopy Several analytical tools have been applied to further the understanding and modeling of hydrogen bonding including Raman scattering [153 ,188 190] , neutron diffraction [191 193] and NMR spectroscopy [13,14,194] . Infrared spectro scopy [146,150,169,195 198] has been especially instrumental in uncovering the behavior of molecules that hydrogen bond. There is an established understanding that the effects of hydrogen bonding are most ev ident in the changes occurring in the O - H vibration region. For this reason, research is often focused on investigating the effects of various conditions, such as temperature and pressure , on the characteristics of the hydroxyl vibrational band [146] . As hydrogen bonded cluster sizes increase, the high frequency O - H vibrational peak decreases and a second broad O - H peak increases at lower frequencies [144,150,154] . The spectrum is fit ted with curves to represent the experimental data. The contributions of the hydroxyl bonds to the band is known in general terms [199] . However, the curve fitting procedure has historically been conducted without recognizing the importance of variations in extinction coefficient with wavenumber . Deconvolution typi cally uses several 152 symmetrical peaks that are a ssigned to specific bonds or species based on the absorbance dependence on temperature or pressure. As indicated by Barlow [200] , a common challenge with this approach is that the problem is under - constrained; comparable fits can be obtained from several different interpretations depending on the initial guess. Computational simulations are increasingly being leveraged to lend more insight to hydrogen bonding in alcohols [17,139,154,155,160] and its effect on their IR vibrations. Studies have included molecular dynamics (MD) [153] , quantum mechanical (QM) [139,144,154 159] and hybrid multi - level (QM/M D ) simulations [17,160 162] .While MD simulations are less computationally expensive than QM calculations, they neglect quantum effects which are important in hydrogen bonding. Another sh ortcoming of MD was highlighted by Kwac and Geva [17] . Using simulations of CH 3 OD in carbon tetrachloride, they found that the relative populations of hydrogen bonded clusters calculated are dependent on the harmonicity, polarizability and damping of the empirical MD forcefield . The relative populations of the different bonds varied by up to 20% percent depending on the choice of forcefield. Furthermore, in their MD simulations with ethanol + water mixtures, Gereben et al. [201,202] were not able to identify a forcefield model that represented experimental scattering X - ray structure factors well across the composition range. Instead, different forcefields were recommended depending on the composition simulated. In Chapter 6 , we developed a sequential MD + QM method to calculate the IR vibrational frequencies and integrated absorption coefficients of different types of O - H bonds . The se wer e classified according to how their atoms participate in hydrogen bonding a s shown in Figure 7 . 1 . 153 Figure 7 . 1 : Types of covalent OH bond s Though all the bonds in Figure 7 . 1 occur in solution, it is often assumed, especially in chemical engineering thermodynamics, that and bonds are formed in minimal quantities and are thus ignored. Moreover, cyclic species are usually neglected. Gereben et al . recently extended their work with ethanol + wate r to study cluster formations [203] and found that 60% (for an 80 mol% ethanol in water mixture) to 95% (pure water) of hydrogen bonded molecules were bonded in cyclic clusters. Considering this, the assumption of no cyclic species is diffic ult to justify. However, the thermodynamic models which use this assumption are often capable of representing experimental data reasonably well [14,68,173] . Therefore, as a first approximation, we ignore cyclic species in this work and defer their analysis to future work. The work covered in Chapter 6 aimed to quantify concentration of alcohol mole cules in an alcohol + alkane mixture using the infrared O - H vibration. To this end, QM calculations were conducted to find relationships between several spectroscopic properties including the vibrational frequency of an O - H bond and NMR shift of its proton . A key finding was the correlation shown in Figure Accepts: 0 Donates: 0 Accepts: 0 Accepts: 1 Accepts: 1 Donates: 1 Donates: 1 Donates: 0 Accepts: 2 Donates: 1 Acce pts: 2 Donates: 0 154 6 . 6 (and repeated in Figure 7 . 2 ) between the integrated absorption coefficient, , and the wavenumber of O - H vibrations. Figure 7 . 2 : Hydroxyl stretching frequencies and integrated absorption coefficients for n - butanol + and ethanol + cyclohexane mixtures calculated from QM simulation s . Diamond, triangle, square and circle markers denote , , , and bond vibrations respectively. The solid line is a line of best fit determined visually It was shown that the relationship is very well modeled by an inverted Gaussian function give n by: 7 . 1 When used to analyze the experimental spectra, the minimum of the Gaussian function, where and bonds fall, was shifted to match the free - end peak wavenumber in the spectrum by adjusting parameter to 3645 cm - 1 . The quantity was used with the Beer - Lambert law to calculate a - H concentration. In this work, we extend this method to calculate the concentrations of and bonds . Further, this cluster distribution 0 200 400 600 800 1000 1200 1400 3200 3300 3400 3500 3600 3700 Integrated Absorption Coefficient (km/mol) Wavenumber (cm - 1 ) 155 is used to calculate association (or hydrogen bonding) parameters for thermodynamic modeling of the alcohol + alkane systems. 7.1.2. Thermodynamic modeling of hydrogen bonding Capturing the effects of hydrogen bonding is important for the thermodynamic modeling of many systems encountered in the chemical industry. In our previous work [96] , we developed a model that calculates deviations of the real solution from an ideal solution reference straight via an activity coefficient term, , which is unity for an ideal solution and non - unity otherwise. Deviations due to hydroge [47 50] , a statistical mechanics appro ach which assigns sites on molecules which act as hydrogen bond oxygen is designated as an acceptor site and the hydrogen as the donor site. Mathematically, this beh avior is captured through: where and are the fraction of sites which remain free (not bonded) in the mixture and when the site host is pure , respectively . Variable is the apparent mole fraction of th e molecule that hosts site B i . T he number of identical sites on a host is given by . For example, the two hydrogen atoms in water are modeled as identical sites due to the symmetry of their host molecule and the assumption that they wil l interact in the same way with their microenvironments, independent of bonding order. Careful consideration on the definition of the term 7 . 3 ), reveals a physical significance rel evant to IR spectroscopy. 7 . 2 156 7 . 3 Here, and are the moles of sites with free sites and the moles of total sites respectively. For alcohols and alcohol + inert (non - hydrogen bonding component) mixtures, represents two sites simultaneously, the fraction of free acceptor sites (oxygen atoms), , and the fraction of f ree donor sites (hydrogen atoms), , at equilibrium. Because each hydrogen bond involves one acceptor and one donor site, the material balance when the alcohol is the only associating species yields . Coupling this with Eq. ( 7 . 3 ) yields: where is the concentration of and bonds and is the apparent concentration of the an association parameter, . For the case considered here where the alcohol is dissolved in a non - associating solvent: where is the mole fraction of alcohol in the mixture and is the solution molar density. Therefore, by using Eq. ( 7 . 4 ) and Eq . ( 1 . 4 ), it is possible to regress spectroscopic measurements of to find the value of the association parameter . The generalized model and form of Eq. ( 3 . 7 ) are provided in our previous w ork [96] . including the cubic - plus - association (CPA) [11,52] and the statistical associating fluid the ory (SAFT) family of equations [9,54 56,95] . In these models, the association parameter is normally found by regressing macroscopic experimental data such as phase equilibria data, vapor pressures 7 . 4 7 . 5 157 and/or saturated liquid densities. However, as shown in Figure 5 . 37 in Chapter 5, t hese methods can yield values of that differ by factors of up to four to five depending on the data used for the fitting. This discrepancy is greater than the expected dependence on solvent . An alternative approach is to derive associa tion parameter values from data collected using fundamental tools such as IR and NMR spectroscopy . Studies that leverage these tools often calculate a concentration - based equilibrium constant, (or a directly related form of the equilibrium constant ) instead of . In these cases, association is modeled using chemical theory which is a [58,96] . 7.1.3. Motivation In this work, we extend the quantitative MD + QM analysis methods developed in Chapter 6 to calculate the concentrations of different types of O - H bonds in solution. Further, the aim is to describe the methodology and provide preliminary results for the calculation of association parameters for the thermodynamic activity coefficient model in Eq. ( 3 . 7 ). This primary method will inv olve interpretation of the infrared spectra of n - butanol + cyclohexane mixtures. Nuclear magnetic resonance (NMR) spectroscopy is used as a complementary tool to confirm findings from IR. 7.2. Experimental m ethods The complete methodology for the computational and IR studies is identical to that in Chapter 6 and details can be found in Sections 6.2 and 6.3 . For the NMR experiments, two protocols were [98] NMR procedure closely and collected data at T=34 C. Samples were prepared volumetrically in a glove bag under a continuous flow of nitrogen gas. Anhydrous n - butanol, cyclohexane and deuterated cyclohexane at purities of 99.8, 99.5% and 99% respectively, were purchased from Sigma - Aldrich and used with no further purification. For 158 dilute samples with an alcohol mole fraction less than 0.03, the solvent was completely replaced with deuterated cyclohexane (Sigma Aldrich isotopic 99.6 atom % D ) for the purpose of diminishing the solvent proton shift s such that no significant overlap occurs with the alcohol O - H proton shift. For samples at compositions greater than 3 mol% n - butanol, a coaxial tube system was used in which the inner tube contained the prepared sample and the outer tube contained deuter ated cyclohexane (for locking) and TMS (as an internal reference). For dilute samples, the outer tube was left empty and the inner tube contained the sample with 0.5 mole% TMS. Data was obtained using a Varian 600 MHz superconducting NMR Spectrometer and t he temperature was calibrated with ethylene glycol. For data collected at other temperatures, a very similar procedure was used but the coaxial tube system was replaced with a simpler single tube configuration and all samples were prepared using deuterate d cyclohexane as the solvent. In these cases, d ata was obtained using a Varian 500 MHz NMR Spectrometer. 7.3. Results and d iscussions In IR spectroscopy, the consequences of hydrogen bonding are most evident in the region of the O - H stretching vibration, in whi ch the bands combine and appear as two peaks between ~ 3 200 and 37 00 cm - 1 . As mentioned in Chapter 6 , the raw spectra are processed to remove solvent bands and the O - H bands are determined to begin and end at local minima on either side of the two peaks. Fu rthermore, the method developed earlier to create more physically significant spectra by scaling the absorbance with the integrated absorption coefficient function in Figure 7 . 2 is repeated here . However, in this c hapter, the integration ranges determined for each individual spectrum by finding the minima on either end of the hydroxyl band. These ranges are provided in Table 7 . 1 . Table 7 . 1 : Frequency range defined as hydroxyl band for each spectrum 159 Temperature ( C) 30 40 50 60 70 0.0469 3055 - 3750 3064 - 375 2 3091 - 3750 3092 - 3750 3153 - 3752 0.0678 3048 - 375 1 3053 - 3751 306 7 - 3751 307 6 375 2 3089 - 375 2 0.0817 30 50 - 375 1 3058 - 3751 306 6 - 3751 307 6 - 3751 3088 - 3751 0.10 0 3049 - 375 2 3050 - 3755 305 6 - 3752 3066 - 3752 3078 - 3752 Although the integrated ranges are not consistent across all temperatures and compositions, the difference between the integrated areas calculated using this method compared to that in Chapter 6 (in which a fixed range was integrated for all spectra) is only 1.2%. In the future, we plan to conduct the analysis described in this chapter over a fixed wavenumber range for a wider variety of systems and conditions. However, the minor change in the integration limits would not affect the preliminary results provided here. Figures 7 . 3 - 7 . 6 show the raw and scaled spectra for n - butanol + cyclohexane at four concentratio ns . There are some features of the O - H vibrational band that are well - documented which are evident in the raw spectra. First, the band consists of two peaks: a sharp, high - frequency peak and a broad, lower - frequency peak. The former peak is often attribute d to free end group protons, and bonds [146,150,153,154,204] , though there is some disagreement on whether contributes significantly to this peak [185,186] . As the temperature is increased, the sharp free end peak increases in absorbance while the broad O - H peak diminishes. This is evidence of shifting bond distributions as hydr ogen bonds break at higher temperature [144,146] and more monomers form in the solution. Moreover, as the proton in a bond becomes more tethered to other s by hydrogen bonding, its vibrations occur at lower frequencies and higher intensities. This is evident in our own work in Figure 7 . 2 as well as that of others [153,205,206] . 160 Figure 7 . 3 : Raw (top) and scaled (bottom) O - H IR band for a 0.0469 mole fraction of n - butanol in cyclohexane . Data was collected at five temperatures between 30 C and 70 C at a pathlength of 0.01092 cm. 161 Figure 7 . 4 : Raw (top) and scaled (bottom) O - H IR band for a 0.0678 mole fraction of n - butanol in cyclohexane . Data was collected at five temperatures between 3 0 C and 70 C at a pathlength of 0.01092 cm 162 Figure 7 . 5 : Raw (top) and scaled (bottom) O - H IR band for a 0. 0817 mole fraction of n - butanol in cyclohexane . Data was collected at five temperatures between 3 0 C and 70 C at a pathlength of 0.01092 cm. 163 Figure 7 . 6 : Raw (top) and scaled (bottom) O - H IR band for a 0.100 mole fraction of n - butanol in cyclohexane . Data was collected at five temperatures between 3 0 C and 70 C at a pathlength of 0.01092 cm. 164 7.3.1. IR band assignments It has been recognized that IR vibrational bands are more accurately attributed to bonds [147,150,172,207] , as first proposed by Hall and Wood [170] , rather than to entire hydrogen bonded clusters. For this reason, we begin peak assignments by studying how the , , and bonds vibrate in the QM simulations. The number of the different species and bond types analyzed to generate the QM calculated distributions in this work is given in Table 7 . 2 . Table 7 . 2 : Number of each species and bond analyzed with QM calculations for the distribution plot in Figure 7 . 7 and Figure 7 . 8 . All species la rger than a monomer have one and one bond. Trimers and tetramers also possess one and two bonds respectively. Ethanol + Cyclohexane n - But anol + Cyclohexane = 0.5 0.1 and 0.5 Species Monomer 859 1019 Dimer 413 383 Trimer 236 105 Tetrame r 223 70 Bonds 859 1019 872 55 8 682 2 45 872 55 8 Figure 7 . 7 shows the distribution of vibrational frequencies of , , and bonds calculated through QM simulations of ethanol + and n - buta nol + cyclohexane mixtures. Findings from the two systems are plotted together because separate plots showed no discernable difference between them. To create the figure, the vibrations of each bond type were binned in 5 cm - 1 wide bins and normalized, and the resulting lines were smoothed with a 7 - bin moving average and normalized to display the relative populations for each bond type as a function of wavenumber. These processing steps were carried out to refine the shape of the distributions and display re lative information about 165 them more clearly. It should be emphasized that only semi - quantitative information should be obtained from this figure from the peak shapes. The peak heights should not be compared between bond types because the data are insufficie nt to reflect the distribution of bond types at the simulated concentrations. Instead, the value of Figure 7 . 7 lies in its ability to reveal approximate peak shapes and a qualitative description of the vibrational frequencies of different bond types relative to one another. In agreement with other studies [146,150,153,154,204] , Figure 7 . 7 shows that and bonds vibrate at the same frequency and show significant overlap . This confirms that the bond contribution to the sharp high - frequency O - H peak should not be ignored. Moreover, the QM calculations show that some of and hydroxyls are foun d at frequencies red - shifted from the high frequency peak. This finding is noteworthy as most previous work has considered only the high frequency peak to represent all the and hydroxyls [150,153] . Figure 7 . 7 : Smoothed normalized distributions of bond types from QM simulations of n - butanol + and ethanol + cyclohexane mixtures . Vibrational frequencies of (n=1878) , (n=1430) , (n=1430) , and (n=927) bonds are binned in 5 wavenumber bins and smoothed with a 7 - bin moving average 166 Figure 7 . 8 : Smoothed normalized distributions of bond types from QM simulations with from dimers separated from others . Vibrationa l frequencies of (n=1878) , (n=1430) , (n=796) , (n=634) and (n=927) bonds are binned in 5 wavenumber bins and smoothed with a 7 - bin moving average Another interesting observation is made when the vibrations of the bond on dimer clusters, is separated from that of bonds on larger oligomers, . As shown in Figure 7 . 8 , vibrations are blue - shifted by about 50 cm - 1 compared to . Inte restingly, bonds on dimers display no similar shift compared to those on larger clusters. We hypothesize that the reason for the shift in is the small size of the dimer molecule. This causes the frequencies of protons in the bonds to be h igher as they have characteristics that are a hybrid of those of the hydrogen it is attached to and those of bonds on larger clusters, . For species larger than dimers, the proton is separated from proton by one or more bonds, so the vibrational mixing lowers its frequency. To our knowledge, this is the first time this distinction has been made and we use this finding to guide the peak fitting and assignment procedure described in the next section. 167 7.3.2. IR peak fitting Quan titative fitting of infrared spectroscopy requires the overlap between constituent peaks to be accurately determined [208] . As alcohol + inert mixtures become more concentrated in alcohol, extensive hydrogen bonding causes the broad O - H peak to have more and more overlap with the free end p eak. Therefore, as indicated by Wandschneider et al . [155] , it becomes difficult to calculate the concentration of free ends when the mole fraction of alcohols exceeds 0.3. For this reason, we limit our study to the mixtures shown in the spe ctra in Figures 7 . 3 - 7 . 6 where the alcohol mole fractions are below 0.100. Based on the line shapes found from the QM c alculated vibrational frequenc y distributions, we fit each scaled experimental IR spectrum with 6 Gaussian peaks: one for , (visually determined to require 2 peaks), , and . Each Gaussian peak, is given by the function: 7 . 6 where , and are fitted parameters. To focus the search for correct parameters, the regression is constrained such that the concentration of bonds is equal to that of bonds. This physical constraint holds as long as no bra nched chains exist in solution as is assumed in our work. The area under a Gaussian curve is given by , 7 . 7 7 . 8 Initial guesses for the wavenumbers of the peak apexes were deter mined from the shifts of the , , and peaks in Figure 7 . 8 r elative to the peak. The upper and lower bounds set for 168 the wavenumbers was 80 cm - 1 from the initial guesses. Figures 7 . 9 - 7 . 12 show the results of the fitting to the spectra at the four lower temperatures for each of three concentrations studied here. Figure 7 . 9 : Peaks fitted to scaled IR spectra for n - butanol + cyclohexane mixture at =0.0 469 . Figures a, b, c, d and e indicate data at T= 30, 40, 50, 60, 70 C respectively 169 Figure 7 . 10 : Peaks fitted to scaled IR spectra for n - butanol + cyclohexane mixture at =0.0 678 . Figures a, b, c, d and e indicate data at T= 30, 40, 50, 60, 70 C respectively 170 Figure 7 . 11 : Peaks fitted to scaled IR spectra for n - butanol + cyclohexane mixture at =0.0 817 . Figures a, b, c, d and e indicate data at T= 30, 40, 50, 60, 70 C respectively 171 Figure 7 . 12 : Peaks fitted to scaled IR spectra for n - butanol + cyclohexane mixture at =0. 100 . Figures a, b, c, d and e in dicate data at T= 30, 40, 50, 60, 70 C respectively 172 The concentration of each of the bond types is calculated by multiplying the areas under the corresponding peaks by the constant relating scaled peak areas to apparent concentrations given in Figure 7 . 13 . The scaling factor is 1/37.47 = 0.0290. Figure 7 . 13 : Scaled band area as a function of apparent n - butanol concentration . The equation for the line of best fit p rovides a convenient method of relating the calculated apparent concentrations from the scaled band area. The hydroxyl species distributions calculated from the scaling factor are given in shown in Figure 7 . 14 . In plots were a gray bar appears, the fitted spectrum was not able to fit the experimental the bars exceed unity indicate that the fitted peaks over predict the prepared concentration value. These slight discrepancies are especially prominent in the plots for the two most dilute sample and are likely due to experimental error in preparing the samples. The exception to this, though, is the fit of the data at T=70 C and which has an unassigned percentage of 13.1%. The peak shape for the raw spectra in Figure 7 . 3 has unexpected features that are indicative of a problem with the experiment, such as evapor ation. y = 37.47x R² = 0.973 10 15 20 25 30 35 40 0.35 0.45 0.55 0.65 0.75 0.85 0.95 Scaled IR Band Area Apparent Concentration (mol/L) 173 Figure 7 . 14 : Bond distributions calculated from IR peak fitting . Colors show ( ) , ( ) , ( ) , ( ), ( ) bond concentrations and unassigned peak area ( ). Fro m top to bottom, figures pertain to samples at 0.0469, 0.0678, 0.0817 and 0.100 mol fraction of n - butanol. 0 0.2 0.4 0.6 0.8 1 Fraction of total alcohol concentration 0 0.2 0.4 0.6 0.8 1 Fraction of total alcohol concentration 0 0.2 0.4 0.6 0.8 1 Fraction of total alcohol concentration 0 0.2 0.4 0.6 0.8 1 Fraction of total alcohol concentration T ( C) 174 In general, however, the plots show several encouraging trends. First, the percent of bonds increases with temperature by an average of 16.1% from 30 C to 70 C. This is consistent with the fact that more hydrogen bonds break as temperature is increased, converting large clusters to smaller oligomers and monomers. For the same reason, the fractions of bonds should and, for the most part, do decr ease with increasing temperatures. The trends for and are less predictable because their concentrations depend on the rate at which larger clusters break down to form smaller clusters (which increases the number of and bonds) compared to the rat e of dimers breaking to form monomers (which decreases the number of and bonds). However, from Figure 7 . 14 , it is observed that the concentrations decrease slightly but consistently for the three most di lute samples, suggesting that the rate of large cluster disintegration to smaller (non - monomer) species is greater. For the highest concentration, the concentrations decrease until the T = 50 C point then increases again. From this, it may be deduce d that the two rates balance at approximately 50 C at this concentration. An interesting trend is noted when the concentration of bonds is plotted as a function of n - butanol mole fraction (see Figure 7 . 15 ): Figure 7 . 15 : Concentration of bonds vs. mole fraction of n - butanol . 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 Concentration of bonds Mole fraction of n - butanol 175 There is a linear decrease at all temperatures except for the data point at T=70 C and . This is further evidence that the data point is flawed and should be measured again and refitted. However, the relationship shown in the figure for the other data points is striking. To calculate the association parameters in the following section, the concentration of free ends must be known. Therefore, the calculated values for are given in Table 7 . 3 . Table 7 . 3 : Concentration of free ends ( ) in mol/mL calculated from fitted Gaussian peaks x 30 40 50 60 70 0.0469 2.1 6E - 04 2.36E - 04 2.51E - 04 2.63E - 04 2.43E - 04 0.0678 3.25E - 04 3.51E - 04 3.70E - 04 3.92E - 04 4.03E - 04 0.0817 3.65E - 04 3.95E - 04 4.14E - 04 4.31E - 04 4.50E - 04 0.1001 4.04E - 04 4.37E - 04 4.61E - 04 4.96E - 04 5.20E - 04 7.3.3. IR calculation of association parameters To extract q uantitative information from IR spectroscopy, it is necessary to assume some form of model for the species in solution and assign bands accordingly. Recently, Chen et al. [209] calculated the dimerization equilibrium constant for 3 - ethyl - 2 - methyl - 3 - pentanol in tetrachloroethylene for quantitative analysis of infrared spectra. This sy stem was intentionally chosen for the study because steric hindrance was thought to limit the extent of association to dimers. The analysis was consistent with this notion, but simple computer simulations point to little trouble going from dimers to linear and cyclic trimers. Such limited hydrogen bonding systems are, of course, unusual among alcohols which hydrogen bond extensively. For association of alcohols, models and band assignments vary widely in literature. A common model for alcohol + inert system s assumes that chains of infinite sizes can form continuously in solution. Such models have been used with some success to interpret associating systems [14,37,68,89,210 212] . 176 In contrast, Asprion et al. [7,137,138] published several papers which assume only monomer, dimer and one other n - mer formation (1 - 2 - n ) where n was a fitted parameter. The work supplement ed an existing thermodynamic model, UNIQU AC, with chemical theory and derived a value of from infrared spectra. Many other studies obtain best fits of their experimental data by assuming at least one cyclic species (most often a cyclic trimer [150,213] or tetramer [144,185,186,214] ) is stable enough to exist in appreciable quantities at equilibrium. In this work, we aim to derive a value fo r the Wertheim association constant, and therefore consider only linear species. First, we regress the values of to best fit t he calculated free end concentrations in Table 7 . 3 using Eq. ( 1 . 4 ). The resulting fit is given in Figure 7 . 16 . Figure 7 . 16 : Calculated and fitted concentrations of free end ( hydroxyls . Markers indica te the concentrations derived from IR at T= 30 ( crosses ), 40 (triangles), 50 (circles), 60 (filled circles), 70 (diamonds) C respectively. Then, the value is further separated using the form used in CPA: 7 . 9 177 The and are regressed to best fit the found in the first step. The quantity in Eq. ( 7 . 9 ) is the radial distribution func tion defined as: 7 . 10 where is the CPA covolume term and is the mixture apparent molar density. The final fit of is giv en in Figure 7 . 17 . Figure 7 . 17 : Derived association parameter values from IR analysis (solid line) compared to CPA values (dashed line) for the self - assoc iation of n - butanol . Markers indicate the values from fitting the IR measurements. While we plot the CPA value of in the figure as well, it should be recognized that the CPA calculations are not the standard to which other values of must be compared. However, it does provide reference values. It can be seen that the association parameter calculated from the IR fits has a weaker temperature dependence with its parameter being 18% smaller than that of CPA. This paramet er is a measure of the enthalpy of the hydrogen bond and the value calculated 178 from IR (2065 K) corresponds to an energy value of 4.1 kcal/mol which is in the correct range for hydrogen bonds. The IR - derived parameter is about 12% larger than that of CPA though the physical significance of is less defined. In the next section, we leverage another spectroscopic tool, NMR, for the calculation of association parameters. 7.3.4. NMR calculations and results A useful complement to IR is proton nuclea r magnetic resonance ( 1 H NMR) spectroscopy. Chemical shifts of proton spins observed in NMR spectra are a result of differences in the electron shielding of hydrogen atoms on a molecule due to their immediate environment. Unlike IR and NIR, the times cale o f an NMR measurement is long relative to the lifetime of a hydrogen bond. Therefore, NMR spectra illustrate a time - averaged observation rather than an instantaneous one. However, there are no issues with varying extinction coefficients as there are with IR . Furthermore, with modern equipment, this form of spectroscopy is highly sensitive [98] and accurate to 0.1 ppb [215] . For alcohol + hydrocarbon systems, the work of Gutowsky and Saike [216] is heavily cited for the derivation of the relationship between the observed shift, , for an O - H site and the weighted average of the shift of the free OH groups and the associated OH groups given by: where , and are the observed chemical shift, the monomer chemical shift and the chemical shift of hydrogen bonded protons respectively. is the concent ration of linear clusters of size i . I t is assumed that end - groups and monomers have the same shift and all complexed species (linear and cyclic produce the same shift. 7 . 11 179 In this section, we derive the values of the self - association parameter for n - butanol from NMR. To this end, we test out two models: the 1 - K Continuous Linear Association Model (CLAM) that was implemented by several studies [68,98,173] and a new model that we have developed, called the 2 - K CLAM . As indicated by its name, the 1 - K CLAM model neglects non - linear species and assu mes that the formation of every hydrogen bond is governed by the same equilibrium constant, . While the simplicity of this type of model is appealing, its often insufficient in capturing real behavior [98] . For this reason, Karachewski et al. developed several more sophisticated models including the LACT [98] (linear association with cyclic trimer) model and AVEC [13] (association with variable equilibrium constant) models. These models require more fitted parameters and, in the case of the AVEC model, have no closed form. This While this is unlikely to affect the physical significance of the model parameters significantly, any a priori assumptions that must be made are a disadvantage. For the 2 - K CLAM model, we recognize an important phenomenon regarding the energetic and entropic effects that accompany the formation of a dimer. The formation of the first hydrogen bond releases less energy than su [155] . When a hydroxyl group is already accepting or donating in a hydrogen bond, the charge redistribution is such that the bond be comes more receptive to another association. As such, it seems justified to include a different for dimer formation than other oligomers. The 2 - K CLAM does this by: 1. Combining a complete mole balance of alcohol molecules with the law of mass action using the equilibrium constant for larger oligomer formation, 2. S ubtracting the dimer term from the mole balance 3. Re - introducing the dimer ter m with a different equilibrium constant , 180 Similar to our calculations with IR spectra, we derive the association parameters from NMR spectra by fitting the CPA equation in Eq. ( 7 . 6 ). Therefore, the 1 - K CLAM requires fitting 3 parameters: , and the chemical shift of hydr ogen - bonded protons, . The monomer shift is obtained experimentally as the shift visible for an extremely dilute n - butanol + cyclohexane mixture ( <0.002). The 2 - K CLAM model requires fitting 5 parameters: , , , and . Figures 7 . 18 - 7 . 20 show the observed chemical shifts relative to the monomer shift for various m ole fractions of n - butanol in cyclohexane at three temperatures. The black lines show the fits obtained by the 1 - K and 2 - K CLAM models when all the parameters, including the parameters, are fitted. Figure 7 . 18 : O - H proton chemical shift as a function of n - butanol mole fraction for n - butanol + cyclohexane mixtures at 26.3 C . Inset shows a close - up of the dilute region. The markers, dashed line and solid line indicate experimental data, the 1 - K CLAM model and the 2 - K CLAM model respectively. The black lines include the parameter in the regression while red lines fix at the CPA value of 2526 K. 181 Figure 7 . 19 : O - H proton chemical s hift as a function of n - butanol mole fraction for n - butanol + cyclohexane mixtures at 34 C . Inset shows a close - up of the dilute region. The markers, dashed line and solid line indicate experimental data, the 1 - K CLAM model and the 2 - K CLAM model respect ively. The black lines include the parameter in the regression while red lines fix at the CPA value of 2526 K. Figure 7 . 20 : O - H proton chemical shift as a function of n - butanol mole fra ction for n - butanol + cyclohexane mixtures at 41 C . Inset shows a close - up of the dilute region. The markers, dashed line and solid line indicate experimental data, the 1 - K CLAM model and the 2 - K CLAM model respectively. The black lines include the parameter in the regression while red lines fix at the CPA value of 2526 K. 182 The shortcomings of the 1 - K CLAM are very evident in both the concentrated and dilute region. In contrast, the 2 - K CLAM is much more capable of describing the dilute region where dimers are a dominant species. The regressed parameters calculated from NMR are given in Table 7 . 4 . Table 7 . 4 : Regressed parameters from NMR analysis and corresponding CPA parameters for comparison CPA IR derived 1 - K CLAM (NMR) 2 - K CLAM (NMR) , , 6.535E - 01 7.31E - 01 3.135 E - 03 1.722 E+00 5.118 E - 02 2.526E+03 2.07E+03 4. 182 E+03 1.712 E+03 3.5 2 5E+03 6.1 36 5.28 1 The parameters are derived from NMR are very different from both CPA values and those derived from the IR fitting. The parameter of the 1 - K CLAM model is almost double that of derived from the IR fits (a difference in the enthalpy of the h ydrogen bond of approximately 4.2 kcal/mol) and the parameter is less than three orders of magnitude smaller . With the 2 - K CLAM, the dimer enthalpy parameter is very low (3.4 kcal/mol) and that of larger oligomers is very high (7 kcal/mol). The N MR algorithm can yield multiple solutions with similar fits to the data. To demonstrate this, we fixed the value of to the CPA value of 2526 K for both the 1 - K CLAM and 2 - K CLAM and refitted the other parameters. This yields the red lines in Figures 7 . 18 - 7 . 20 . The values of the shift remained almost the same but , and were found to be 0.696, 0.1166 and 1.341 res pectively. These quantities are closer to CPA and IR - derived values. This indicates that a more reasonable solution exists but that the regression program must be constrained in order to find it. More work must be done to correct for this and we outline th is and other future work in Section 7. 4 . 183 7.4. Summary of future work This chapter contains preliminary analyses and results for the calculation of Wertheim association parameters from IR and NMR spectroscopy. As an ongoing project, there are several areas that will be further pursued in the future. Firstly, the fitting of the scaled infrared spectra should be revisited. Specifically, the spectrum at T=70 C and produced strange trends as discu ssed earlier. This may be due to errors in the experimental IR data collection because evaporation can occur at high temperatures. However, this may also be due to the low wavenumber cutoff which marks the end of the integration area. This is most clearly visible in Figure 7 . 3 where the line is discontinued at a lower wavenumber than the other spectra. Currently, the range of the plot is determined by finding the minima on either side of the O - H band. However, in th e future, data should be integrated between band limits that are consistent for all temperatures and concentrations. Though the QM species distributions in Figure 7 . 8 su ggest that the hydroxyl, like the hydr oxyl, is best described by two Gaussian peaks instead of one, adding another peak to the fitting procedure may result in difficulties in discerning differences between the two types. For this work, the hydroxyls were represented by a broad low peak in the high frequency region and the hydroxyls were represented by a tall narrow peak. However, because the and distributions are so similar in Figure 7 . 8 separation of the species may not be correct. Rather, it may be more realistic combine the two peaks together as the sum , eliminate the area constraint between and hydroxyls, and, instead. focus on improved initial guesses for the hydroxyls. Regarding the shape of the fitted peaks, it may be w orthwhile to e xplore other line shapes beyond Gaussian such as Lorentzian or gamma distribution functions. Some researchers have found success by combining two functions such as with the Gaussian - Lorentzian curves [200] . However, 184 this requires including more parameters for unclear benefit. Instead, it is recommended that Lorentzian functi ons be used for the peaks that convolute to form the sharp high frequency peak and Gaussians be used for the others. In this way, the same total number of parameters is maintained but the two peaks are analyzed differently. Secondly, we encountered some c hallenges in calculating association parameters from NMR spectra. The key issue is that many combinations of parameters can yield reasonably accurate representations of the experimental data. In such cases, it is important that the stability of the final s olution be tested by perturbing the initial guesses and ensuring that the same result is reached [200] . To correct for this is sue, it is recommended that more constraints be defined to remove degrees of freedom. Moreover, although the 2 - K CLAM model was shown to represent NMR data more accurately, the accuracy of 2 - K Therefore, this endeavor warrants further study. Lastly, none of the models tested in this work considered cyclic species. With linear species, an n - mer contains n - 1 hydrogen bonds whereas a cyclic n - mer contains only n hydrogen bonds. This counting discrepancy is not accounted for with difficult to adapt the function pursuing because, if successful, the resulting model will be a valuable contribution to the field. Additionally, hydrogen bonding in these types of clusters can be studied using the MD+QM m ethod developed in Chapter 6 . This insight can help improve chemical theory models by providing a physical basis for the assumptions used. 185 APPENDIX 186 Appendix I : Derivation of the 2 - K CLAM model To derive the 2 - K CLAM model, we start with an apparent mole balance written in terms of concentration by dividing both sides by total volume, Defining two equilibrium constants, one for the formation of dimers and another for all subsequent oligomers: Combining Eqs. ( I. 3 and I. 4 ) and rearrang ing yields Eq. ( I. 5 ) can be substituted into the mole balance in Eq. ( I. 2 ), Rearranging: Recognizing the converging series I. 8 Eq. ( I. 7 ) can be written as I. 1 I. 2 I. 3 I. 4 for I. 5 I. 6 I. 7 187 Substituting and recognizing that the apparent molar volume of an alcohol, is the volume of the monomer yields Rearranging in terms of , Eq. ( I. 12 ) can be written in the form of a cubic equation where , , and are: The chemical shift is then calculated from Eq. (7.11): Substituting Eq. ( I. 2 ) for and rearranging yields: I. 9 I. 10 I. 11 I. 12 I. 13 I. 14 I. 15 I. 16 I. 17 I. 18 188 or where and represent the observed chemical shift and the chemical shift of hydrogen - bonded protons relative to the monomer shift respectively. Written in terms of volume fractions and assuming that t he volume of an - mer is given by : Finally, substituting the defi nitions of the equilibrium constants in terms of volume fractions (Eqs. ( I. 3 and I. 4 )) and rearranging gives: I. 19 I. 20 I. 21 I. 22 189 Chapter 8. Conclusions and Future Directions The objective of this project is to combine spectroscopic tools and computational techniques to improve the thermodynamic mod eling of systems with hydrogen bonding. To this end, this document begins by outlining two common methods for modeling association in thermodynamics: concentratio n - and activity - based equilibrium constants for the liquid, we first explore this question in C theory [30] and that the concentration - based equilibrium constant, , can only be assum ed to be independent of composition when residual contributions to non - ideality are also constant. In C hapter 3, , which is the association parameter in chemical theory, is shown to be identical to , the association parameter in We reasonable assumptions of zero excess volume and universal packing factor. A Wertheim activity coefficient model (WAG) is also introduced based on methods developed by Michelsen et al . [66] . Furthermore, algebraic and numerical proofs of the equivalence of the two theories are derived for these systems when only linear species are present. This finding is extended to other systems that cont ain both associating and solvating components in C hapter 4. The activity coefficient model developed in C hapter 3 is tested in C hapter 5 using literature (CPA) values and its performance is compared to that of existing association equations of state. It is found that, when are able to capture the behavio r of methanol and ethanol - containing systems well, generating fits with the lowest overall SSQ value. Among the two parameter models, WAG - SH, which combines WAG with Scatchard - Hildebrand provided the best fits and predictions of the experimental data. Inde ed, WAG - SH performed similarly or better than the association equations of state tested, PC - 190 SAFT and CPA. Chapters 6 and 7 describe work on the analysis of IR spectra for the calculation of values. In the former, the hydroxyl vibrational region in IR spectra of alcohol + alkane systems is quantified for the first time and related to the apparent concentration of alcohol in the mixture. This is accomplished by using relations hips between IR properties found by MD + QM calculations of hydrogen bonded clusters. In C hapter 7, preliminary calculations for the value of from the analysis of IR spectra of n - butanol + cyclohexane mixtures are shown. This deconvoluti on procedure lies beyond the scope of this thesis. However, detailed recommendations for future directions in the IR analysis portion of this project are outlined. Other future work is categorized and summarized in Table 8 . 1 . Table 8 . 1 : Accomplished and future g oals for research project Project Goal Completed Future Work Thermodynamic Model Explore similarities and differences between association modeling theories Develop a n activity coefficient model for (WAG model) Test modeling capabilities of WAG for VLE and LLE of alcohol + inert systems Test predictive power for LLE of ternary systems of alcohol + inert systems Test modeling cap abilities of WAG for systems with association and solvation Test modeling capabilities of WAG for components that only solvate Develop a 2 - parameter WAG model analogous to 2 - K CLAM Create an AspenPlus® user model for industrial partners to expl ore Release AspenPlus® user model for public use Spectroscopy Establish reliable IR protocols and benchmark methods against literature Establish reliable NMR protocols and benchmark methods against literature Assemble database of spectroscop ic data for alcohol + inert systems Assemble database of spectroscopic data for other challenging systems 191 Table 8 . 1 Project Goal Completed Future Work MD + QM Develop scripts and run MD simulat ions Develop scripts to analyze MD simulations for hydrogen bonding information Develop scripts to conduct QM calculations on linear clusters Develop scripts to conduct QM calculations on non - linear clusters Run MD+QM calculations for alco hol + inert systems Run MD+QM calculations for other challenging systems Overall Use MD+QM findings to quantify IR peak Use MD+QM findings to quantify concentration of free ends Calculate association parameters from spectroscopy (IR, NMR) an d MD + QM Test modeling and predictive power of WAG with new association parameters For the thermodynamic model development, the derivation of the WAG function should be revisited to relax some of its underlying assumptions such as the use of a con stant packing fraction. Moreover, the WAG model should be tested for more systems that contain solvating components. An especially complex case occurs when components that do not associate but are capable of solvating, such as dioxane and chloroform, are m ixed. In these cases, the combining rules introduced in C hapter 5 cannot be used and association parameters must be estimated through other means. Furthermore, for consistency across all phases, it will be beneficial to extend the WAG model to be able to c model [79] . Once this is completed, the WAG user model, which has already been developed to function with AspenPlus®, a process design softwa re package, can be released for commercial use. In this work, we developed a chemical theory variation with two equilibrium constants: one for the formation of a dimer and another for subsequent association. The model is shown to be 192 more capable of fitting theory is more easily generalizable than chemical theory. Therefore, future efforts should be invested in developing a Wertheim association theory that allows for two association parame ters. In terms of computational experimentation, future work should include MD + QM calculations like those developed in C hapter 6 for different classes of systems. Additionally, the characteristics of non - linear clusters, such as ring and branched structu res, should be investigated. 193 BIBLIOGRAPHY 194 BIBLIOGRAPHY [1] G.M. Kontogeorgis, M.L. Michelsen, G.K. Folas, S. Derawi, N. von Solms, E.H. Stenby, Ten Years with the CPA (Cubic - Plus - Association) Equation of State. Part 1. Pure Compounds and Self - Associating Systems, Ind. Eng. Chem. Res. 45 (2006) 4855 4868. doi:10.1021/ie051305v. [2] G.M. Kontogeorgis, I. Tsivintzelis, N. von Solms, A. Grenner, D. Bogh, M. Frost, A. Knage - Rasmuss en, I.G. Economou, Use of monomer fraction data in the parametrization of association theories, Fluid Phase Equilibria. 296 (2010) 219 229. doi:10.1016/j.fluid.2010.05.028. [3] H. Renon, J.M. Prausnitz, Local compositions in thermodynamic excess functions for liquid mixtures, Aiche J. 14 (1968) 135 - . doi:10.1002/aic.690140124. [4] L.S. Budantseva, T.M. Lesteva, M.S. Nemtsov, Liquid - vapor equilibrium in systems methanol - C7 - 8 hydrocarbons of different classes, Zh. Fiz. Khim. 49 (1975) 1844. [5] H. Katayama, L iquid liquid equilibria of two ternary systems: methanol cyclohexane including 1,3 - dioxolane or 1,4 - dioxane in the range of 277.79 308.64 K, Fluid Phase Equilibria. 164 (1999) 83 95. doi:10.1016/S0378 - 3812(99)00241 - 1. [6] Y. - H. Fu, S.I. Sandler, H. Orbey, A Modified UNIQUAC Model That Includes Hydrogen Bonding, Industrial & Engineering Chemistry Research. 34 (1995) 4351 4363. doi:10.1021/ie00039a027. [7] N. Asprion, H. Hasse, G. Maurer, Application of IR - spectroscopy in thermodynamic investigations of assoc iating solutions, Fluid Phase Equilibria. 205 (2003) 195 214. doi:10.1016/S0378 - 3812(02)00181 - 4. [8] J.R. Elliott, C.T. Lira, Introductory Chemical Engineering Thermodynamics, International Edition, 2nd, International ed., Pearson Education, Inc., Upper Sa ddle River, NJ, 2012. [9] W.G. Chapman, K.E. Gubbins, G. Jackson, M. Radosz, New reference equation of state for associating liquids, Industrial & Engineering Chemistry Research. 29 (1990) 1709 1721. doi:10.1021/ie00104a021. [10] J.R. Elliott, S.J. Suresh, M.D. Donohue, A simple equation of state for nonspherical and associating molecules, Ind. Eng. Chem. Res. 29 (1990) 1476 1485. doi:10.1021/ie00103a057. [11] G.M. Kontogeorgis, E.C. Voutsas, I.V. Yakoumis, D.P. Tassios, An Equation of State for Associating Fluids, Industrial & Engineering Chemistry Research. 35 (1996) 4310 4318. doi:10.1021/ie9600203. 195 [12] D.S. Abrams, J.M. Prausnitz, Statisical Thermodynamics of Liquid Mixtures: A New Expression for Excess Gibbs Energy of Partly or Completely Miscible Syst ems, AIChE Journal. 21 (1975) 116 128. doi:10.1002/aic.690210115. [13] A.M. Karachewski, W.J. Howell, C. a. Eckert, Development of the AVEC model for associating mixtures using NMR spectroscopy, AIChE Journal. 37 (1991) 65 73. doi:10.1002/aic.690370106. [1 4] A.M. Karachewski, M.M. McNiel, C.A. Eckert, A study of the hydrogen - bonding in alcohol - solutions using NMR spectroscopy, Industrial & Engineering Chemistry Research. 28 (1989) 315 324. doi:10.1021/ie00087a010. [15] I. Nagata, K. Tamura, K. Gotoh, Associ ation model of fluids. Phase equilibria and excess enthalpies in mixtures containing alcohol and acetonitrile, Z. Phys. Chemie - Int. J. Res. Phys. Chem. Chem. Phys., Frankfort. 199 (1997) 1 23. [16] S.W. Campbell, Chemical theory for mixtures containing any number of alcohols, Fluid Phase Equilibria. 102 (1994) 61 84. [17] K. Kwac, E. Geva, A mixed quantum - classical molecular dynamics study of the hydroxyl stretch in methanol/carbon tetrachloride mixtures: equilibrium hydrogen - bond structure and dynamics at the ground state and the infrared absorption spectrum., The Journal of Physical Chemistry. B. 115 (2011) 9184 94. doi:10.1021/jp204245z. [18] G.M. Wilson, Vapor - Liquid Equilibrium. XI. A New Expression for the Excess Free Energy of Mixing, Journal of the A merican Chemical Society. 86 (1964) 127 130. [19] G. Scatchard, Equilibria in Non - electrolyte Solutions in Relation to the Vapor Pressures and Densities of the Components., Chem. Rev. 8 (1931) 321 333. doi:10.1021/cr60030a010. [20] J.H. Hildebrand, S.E. Wo od, The Derivation of Equations for Regular Solutions, The Journal of Chemical Physics. 1 (1933) 817 822. doi:10.1063/1.1749250. [21] H.S. Fogler, Elements of Chemical Reaction Engineering, Pearson Education, 2016. https://books.google.com/books?id=uQ3rsgE ACAAJ. [22] O. Levenspiel, Chemical Reaction Engineering, Wiley, 1999. https://books.google.com/books?id=X4svAQAAIAAJ. [23] D.A. McQuarrie, J.D. Simon, Physical Chemistry: A Molecular Approach, University Science Books, 1997. [24] K.G. Denbigh, The Princip les of Chemical Equilibrium: With Applications in Chemistry and Chemical Engineering, Cambridge University Press, 1981. https://books.google.com/books?id=cwWT48wL238C. [25] S.I. Sandler, Chemical, Biochemical, and Engineering Thermodynamics, John Wiley & S ons, 2006. https://books.google.com/books?id=4MXDAgAAQBAJ. 196 [26] - action kinetics, Progress in Reaction Kinetics and Mechanism. 30 (2005) 3 113. doi:10.3184/007967405777874868. [27] R. Rönnback, T. Salmi, A. V uori, H. Haario, J. Lehtonen, A. Sundqvist, E. Tirronen, Development of a kinetic model for the esterification of acetic acid with methanol in the presence of a homogeneous acid catalyst, Chemical Engineering Science. 52 (1997) 3369 3381. doi:10.1016/S0009 - 2509(97)00139 - 5. [28] A.K. Kolah, N.S. Asthana, D.T. Vu, C.T. Lira, D.J. Miller, Reaction kinetics of the catalytic esterification of citric acid with ethanol, Industrial and Engineering Chemistry Research. 46 (2007) 3180 3187. doi:10.1021/ie060828f. [29] A. Fredenslund, R.L. Jones, J.M. Prausnitz, Group - contribution estimation of activity coefficients in nonideal liquid mixtures, AIChE Journal. 21 (1975) 1086 1099. doi:10.1002/aic.690210607. [30] P.J. Flory, Thermodynamics of Heterogeneous Polymers and Th eir Solutions, The Journal of Chemical Physics. 12 (1944) 425. doi:10.1063/1.1723887. [31] G.K. Folas, G.M. Kontogeorgis, M.L. Michelsen, E.H. Stenby, Application of the cubic - plus - association equation of state to mixtures with polar chemicals and high pre ssures, Industrial & Engineering Chemistry Research. 45 (2006) 1516 1526. [32] A. Apelblat, The concept of associated solutions in historical development. Part 1. The 1884 1984 period, Journal of Molecular Liquids. 128 (2006) 1 31. doi:10.1016/j.molliq.200 6.02.005. [33] 1. The 1884 - 162. [34] F. Dolezalek, Zur Theorie der Binären Gemische und Konzentrierten Lösungen, Z. Phys. Chem. 64 (1908) 727 747. [35] I. Nagata, K. Tamura, S. Tokuriki, Excess enthalpies and complex formation of acetonitrile with acetone, chloroform and benzene, Thermochimica Acta. 47 (1981) 315 331. [36] I. Nagata, Y. Kawamura, Thermodyna mics of Alcohol - Unassociated Active Component Liquid Mixtures, Chemical Engineering Science. 34 (1979) 601 611. [37] N.D. Coggeshall, E.L. Saier, Infrared Absorption Study of Hydrogen Bonding Equilibria, Journal of the American Chemical Society. 73 (1951) 5414 5418. doi:10.1021/ja01155a118. [38] E.A. Guggenheim, Statistical thermodynamics of mixtures with zero energies of mixing, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences. 183 (1944) 203 212. [39] P. Muthukumara n, R.L. Brinkley, R.B. Gupta, Lattice - fluid equation of state with hydrogen - bond cooperativity, Aiche Journal. 48 (2002) 386 392. doi:10.1002/aic.690480221. 197 [40] R.B. Gupta, R.L. Brinkley, Hydrogen - bond cooperativity in 1 - alkanol plus n - alkane binary mixtures, Aiche Journal. 44 (1998) 207 213. doi:10.1002/aic.690440122. [41] C. Panayiotou, M. Pantoula, E. Stefanis, I. Tsivintzelis, I.G. Economou, Nonrandom Hyd rogen - Bonding Model of Fluids and Their Mixtures. 1. Pure Fluids, Ind. Eng. Chem. Res. 43 (2004) 6592 6606. doi:10.1021/ie040114+. [42] C. Panayiotou, I. Tsivintzelis, I.G. Economou, Nonrandom Hydrogen - Bonding Model of Fluids and Their Mixtures. 2. Multico mponent Mixtures, Ind. Eng. Chem. Res. 46 (2007) 2628 2636. doi:10.1021/ie0612919. [43] K.E. Gubbins, Perturbation Theories of the Thermodynamics of Polar and Associating Liquids: A Historical Perspective, Fluid Phase Equilibria. (2015) 1 15. doi:10.1016/j .fluid.2015.12.043. [44] P.T. Cummings, G. Stell, Statistical mechanical models of chemical reactions, Molecular Physics. 51 (1984) 253 287. doi:10.1080/00268978400100191. [45] G. Stell, Y. Zhou, Chemical association in simple models of molecular and ionic fluids, The Journal of Chemical Physics. 91 (1989) 3618 3623. doi:10.1063/1.456894. [46] Y. Zhou, G. Stell, Chemical association in simple models of molecular and ionic fluids. II. Thermodynamic properties, The Journal of Chemical Physics. 96 (1992) 1504 1506. doi:10.1063/1.462872. [47] M.S. Wertheim, Fluids with highly directional attractive forces.1. Statistical thermodynamics, Journal of Statistical Physics. 35 (1984) 19 34. doi:10.1007/bf01017362. [48] M.S. Wertheim, Fluids with highly directional attr active forces. 2. Thermodynamic perturbation - theory and integral - equations, Journal of Statistical Physics. 35 (1984) 35 47. doi:10.1007/bf01017363. [49] M.S. Wertheim, Fluids with highly directional attractive forces. 3. Multiple attraction sites, Journal of Statistical Physics. 42 (1986) 459 476. doi:10.1007/bf01127721. [50] M.S. Wertheim, Fluids with highly directional attractive forces. 4. Equilibrium polymerization, Journal of Statistical Physics. 42 (1986) 477 492. doi:10.1007/bf01127722. [51] S.J. Su resh, J.R. Elliott, Multiphase Equilibrium Analysis via a Generalized Equation of State for Associating Mixtures, Industrial & Engineering Chemistry Research. 31 (1992) 2783 2794. [52] I. V Yakoumis, G.M. Kontogeorgis, E.C. Voutsas, E.M. Hendriks, D.P. Tas sios, Prediction of Phase Equilibria in Binary Aqueous Systems Containing Alkanes , Cycloalkanes , and Alkenes with the Cubic - plus - Association Equation of State, Industrial & Engineering Chemistry Research. 37 (1998) 4175 4182. 198 [53] J. Gross, G. Sadowski, Perturbed - Chain SAFT: An Equation of State Based on a Perturbation Theory for Chain Molecules, Industrial & Engineering Chemistry Research. 40 (2001) 1244 1260. doi:10.1021/ie0003887. [54] Y. - H. Fu, S.I. Sandler, A Simplified SAFT Equation of State for As sociating Compounds and Mixtures, Industrial & Engineering Chemistry Research. 34 (1995) 1897 1909. doi:10.1021/ie00044a042. [55] S.H. Huang, M. Radosz, Equation of state for small, large, polydisperse, and associating molecules, Industrial & Engineering C hemistry Research. 29 (1990) 2284 2294. doi:10.1021/ie00107a014. [56] A. Gil - Villegas, A. Galindo, P.J. Whitehead, S.J. Mills, G. Jackson, A.N. Burgess, Statistical associating fluid theory for chain molecules with attractive potentials of variable range, The Journal of Chemical Physics. 106 (1997) 4168 4186. doi:10.1063/1.473101. [57] T. Kraska, K.E. Gubbins, Phase Equilibria Calculations with a Modified SAFT Equation of State. 1. Pure Alkanes, Alkanols, and Water, Industrial & Engineering Chemistry Resear ch. 35 (1996) 4727 4737. doi:10.1021/ie9602320. [58] I.G. Economou, M.D. Donohue, Chemical, quasi - chemical and perturbation theories for associating fluids, AIChE Journal. 37 (1991) 1875 1894. doi:10.1002/aic.690371212. [59] J.P. Wolbach, S.I. Sandler, Usi ng molecular orbital calculations to describe the phase behavior of hydrogen - bonding fluids, Industrial & Engineering Chemistry Research. 36 (1997) 4041 4051. doi:10.1021/ie9607255. [60] E.M. Hendriks, J. Walsh, A.R.D. van Bergen, A general approach to ass ociation using cluster partition functions, Journal of Statistical Physics. 87 (1997) 1287 1306. doi:10.1007/BF02181285. [61] J.P. Wolbach, S.I. Sandler, Using molecular orbital calculations to describe the phase behavior of cross - associating mixtures, Ind ustrial & Engineering Chemistry Research. 37 (1998) 2917 2928. doi:10.1021/ie970781l. [62] J.M. Prausnitz, F.W. Tavares, Thermodynamics of Fluid - Phase Equilibria for Standard Chemical Engineering Operations, AIChE Journal. 50 (2004) 739 761. doi:10.1002/ai c.10069. [63] M.L. Michelsen, Robust and efficient solution procedures for association models, Industrial & Engineering Chemistry Research. 45 (2006) 8449 8453. doi:10.1021/ie060029x. [64] G.M. Kontogeorgis, G.K. Folas, Thermodynamic Models for Industrial Applications, John Wiley & Sons, 2010. [65] thermodynamic perturbation theory, Fluid Phase Equilibria. 428 (2016) 121 152. doi:10.1016/j.fluid.2016.07.033. 199 [66] M.L. Michelsen, E. M. Hendriks, Physical properties from association models, Fluid Phase Equilibria. 180 (2001) 165 174. doi:10.1016/S0378 - 3812(01)00344 - 2. [67] J.R. Elliott, C.T. Lira, Errata: Introductory chemical engineering thermodynamics, 2nd ed., 2012, Errata, (2016). http://chethermo.net/errata (accessed July 13, 2016). [68] H. Renon, J.M. Prausnitz, On the thermodynamics of alcohol - hydrocarbon solutions, Chemical Engineering Science. 22 (1967) 299 307. doi:10.1016/0009 - 2509(67)80116 - 7. [69] H. Renon, J.M. Prausnitz, O n the thermodynamics of alcohol - hydrocarbon solutions. Errata., Chemical Engineering Science. 22 (1967) 1891. doi:10.1016/0009 - 2509(67)80223 - 9. [70] I. Nagata, Thermodynamics of alcohol solutions. Phase equilibria and excess molar enthalpies of mixtures co ntaining two alcohols, Thermochimica Acta. 107 (1986) 199 217. [71] A. Nath, E. Bender, On the thermodynamics of associated solutions. I. An analytical method for determining the enthalpy and entropy of association and equilibrium constant for pure liquid substances, Fluid Phase Equilibria. 7 (1981) 275 287. doi:10.1016/0378 - 3812(81)80012 - X. [72] I. Nagata, K. Gotoh, Thermodynamics of alcohol solutions, Thermochimica Acta. 258 (1995) 77 107. [73] V. Brandani, A continuous linear association model for determ ining the enthalpy of hydrogen - bond formation and the equilibrium association constant for pure hydrogen - bonded liquids, Fluid Phase Equilibria. 12 (1983) 87 104. [74] V. Brandani, F. Evangelista, The UNIQUAC associated - solution theory: vapor - liquid equili bria of binary systems containing one associating and one inert or active component, Fluid Phase Equilibria. 17 (1984) 281 302. [75] R. a Heidemann, J.M. Prausnitz, A van der Waals - type equation of state for fluids with associating molecules., Proceedings of the National Academy of Sciences of the United States of America. 73 (1976) 1773 6. [76] A.G. Pradhan, V.R. Bhethanabotla, S.W. Campbell, Vapor - liquid equilibrium data for ethanol - n - heptane - 1 - propanol and ethanol - n - heptane - 2 - propanol and their interpret ation by a simple association model, Fluid Phase Equilibria. 84 (1993) 183 206. doi:10.1016/0378 - 3812(93)85123 - 4. [77] G.M. Lobien, J.M. Prausnitz, Correlation for the ratio of limiting activity coefficients for binary liquid mixtures, Fluid Phase Equilibr ia. 8 (1982) 149 160. doi:10.1016/0378 - 3812(82)80032 - 0. [78] W.A. Fouad, L. Wang, A. Haghmoradi, S.K. Gupta, W.G. Chapman, Understanding the Thermodynamics of Hydrogen Bonding in Alcohol - Containing Mixtures: Self Association, The Journal of Physical Chemis try B. 119 (2015) 14086 14101. doi:10.1021/acs.jpcb.5b08285. 200 [79] Industrial & Engineering Chemistry Process Design and Development. 14 (1975) 209 216. doi:10.1021 /i260055a003. [80] D. Ambrose, J.H. Ellender, C.H.S. Sprake, R. Townsend, Thermodynamic properties of organic oxygen compounds: Vapour pressures of some ethers, Journal of Chemical Thermodynamics. 8 (1976) 165 178. [81] A. Grenner, G.M. Kontogeorgis, N. vo n Solms, M.L. Michelsen, Application of PC - SAFT to glycol containing systems PC - SAFT towards a predictive approach, Fluid Phase Equilibria. 261 (2007) 248 257. doi:10.1016/j.fluid.2007.04.025. [82] J. Suresh, E.J. Beckman, Prediction of liquid - liquid equ ilibria in ternary mixtures from binary data, Fluid Phase Equilibria. 99 (1994) 219 240. doi:10.1016/0378 - 3812(94)80033 - 2. [83] E.C. Voutsas, I.V. Yakoumis, D.P. Tassios, Prediction of phase equilibria in water/alcohol/alkane systems, Fluid Phase Equilibri a. 158 160 (1999) 151 163. [84] I. Nagata, K. Tamura, Thermodynamics of Solutions of Acetonitrile with Propanols, Thermochimica Acta. 98 (1986) 147 158. [85] C.A. Schmidt, S.W. Campbell, An association model for mixtures containing any number of alkanes an d any number of alcohols: Application to alcohol - alkane binary systems, International Journal of Thermophysics. 16 (1995) 1299 1307. doi:10.1007/BF02081297. [86] W. Acree, Thermodynamic properties of nonelectrolyte solutions, Academic Press Inc., 1984. [87 ] J.M. Prausnitz, R.N. Lichtenthaler, E.G. de Azevedo, Molecular thermodynamics of fluid - phase equilibria, 2nd Edition, Pearson Education, 1986. [88] J. Gmehling, D.D. Liu, J.M. Prausnitz, High - pressure vapor liquid equilibria for mixtures containing one o r more polar components: Application of an equation of state which includes dimerization equilibria, Chemical Engineering Science. 34 (1979) 951 958. doi:10.1016/0009 - 2509(79)85006 - X. [89] G.D. Ikonomou, M.D. Donohue, Extension of the associated perturbed anisotropic chain theory to mixtures with more than one associating component, Fluid Phase Equilibria. 39 (1988) 129 159. doi:10.1016/0378 - 3812(88)85002 - 7. [90] I. Nagata, K. Tamura, S. Tokuriki, Excess enthalpies and complex formation of acetonitrile with acetone, chloroform, and benzene, Thermochimica Acta. 47 (1981) 315 331. doi:10.1016/0040 - 6031(81)80110 - 4. [91] W.J. Howell, C.T. Lira, C.A. Eckert, Linear chemical - physical theory model for liquid metal solution thermodynamics, AIChE Journal. 34 (1988) 1 477 1485. doi:10.1002/aic.690340909. 201 [92] I. Nagata, Thermodynamic properties of alcohol - inert component mixtures, Zeitschrift Für Physikalische Chemie, Leipzig. 259 (1978) 1151 21159. [93] I. Nagata, K. Gotoh, Association model of fluids. Phase equilibria in sulphur dioxide mixtures, Fluid Phase Equilibria. 120 (1996) 69 80. [94] I. Nagata, K. Gotoh, Thermodynamics of amine and acetonitrile solutions, Thermochimica Acta. 274 (1996) 81 99. [95] J. Gross, G. Sadow ski, Perturbed - Theory for Chain Molecules, Ind. Eng. Chem. Res. 40 (2001) 1244 1260. doi:10.1021/ie0003887. [96] A.M. Bala, C.T. Lira, Relation of Wertheim association constants to concentration - bas ed equilibrium constants for mixtures with chain - forming components, Fluid Phase Equilibria. 430 (2016) 47 56. doi:10.1016/j.fluid.2016.09.015. [97] D. Ghonasgi, W.G. Chapman, Theory and simulation for associating fluids with four bonding sites, Molecular Physics. 79 (1993) 291 311. doi:10.1080/00268979300101221. [98] A.M. Karachewski, M.M. McNiel, C. a. Eckert, A study of hydrogen bonding in alcohol solutions using NMR spectroscopy, Industrial & Engineering Chemistry Research. 28 (1989) 315 324. doi:10.102 1/ie00087a010. [99] A. Santhanakrishnan, A. Shannon, L. Peereboom, C.T. Lira, D.J. Miller, Kinetics of mixed ethanol/n - butanol esterification of butyric acid with amberlyst 70 and p - toluene sulfonic acid, Industrial and Engineering Chemistry Research. 52 ( 2013) 1845 1853. doi:10.1021/ie302267s. [100] T.L. Jordison, C.T. Lira, D.J. Miller, Condensed - phase ethanol conversion to higher alcohols, Ind. Eng. Chem. Res. 54 (2015) 10991 11000. doi:10.1021/acs.iecr.5b02409. [101] T.L. Jordison, L. Peereboom, D.J. Mi ller, Impact of Water on Condensed Phase Ethanol Guerbet Reactions, Ind. Eng. Chem. Res. 55 (2016) 6579 6585. doi:10.1021/acs.iecr.6b00700. [102] F.E. Anderson, J.M. Prausnitz, Inhibition of gas hydrates by methanol, AIChE Journal. 32 (1986) 1321 1333. doi :10.1002/aic.690320810. [103] H.G. Harris, J.M. Prausnitz, Thermodynamics of Solutions with Physical and Chemical Interactions. Solubility of Acetylene in Organic Solvents, Ind. Eng. Chem. Fund. 8 (1969) 180 188. doi:10.1021/i160030a001. [104] Y. - H. Fu, H. Orbey, S.I. Sandler, Prediction of Vapor - Liquid Equilibria of Associating Mixtures with UNIFAC Models That Include Association, Industrial & Engineering Chemistry Research. 35 (1996) 4656 4666. doi:10.1021/ie950545f. 202 [105] A. Vahid, N.H. Gray, J.R. Elli ott, Trends in the Athermal Entropy of Mixing of Polymer Solutions, Macromolecules. 47 (2014) 1514 1531. doi:10.1021/ma402593m. [106] A.J. Staverman, The entropy of high polymer solutions.Generalization of formulae, Recueil Des Travaux Chimiques Des Pays - B as. 69 (1950) 163 174. [107] M.D. Donohue, J.M. Prausnitz, Combinatorial Entropy of Mixing Molecules that Differ in Size and Shape. A Simple Approximation for Binary and Multicomponent Mixtures, Can. J. Chem. 53 (1975) 1586 1592. doi:10.1139/v75 - 224. [108] I. Kikic, P. Alessi, P. Rasmussen, A. Fredenslund, On the combinatorial part of the UNIFAC and UNIQUAC models, The Canadian Journal of Chemical Engineering. 58 (1980) 253 258. doi:10.1002/cjce.5450580218. [109] U. Weidlich, J. Gmehling, A modified UNIFAC model. 1. Prediction of VLE, hE, and gamma infinity, Industrial & Engineering Chemistry Research. 26 (2012) 1372 1381. doi:10.1021/ie00067a018. [110] U. Weidlich, Experimentelle und theoretische Untersuchungen zur Erweiterung der Gruppenbeitragsmethode UNI FAC, (1985). [111] B.L. Larsen, P. Rasmussen, A. Fredenslund, A Modified UNIFAC Group - Contribution Model for Prediction of Phase Equilibria and Heats of Mixing, Industrial and Engineering Chemistry Research. 26 (1987) 2274 2286. [112] B.L. Larsen, Predicti ons of phase equilibria and heat effects of mixing with a modified UNIFAC model, PhD Thesis, Danmarks tekniske Højskole, 1986. [113] S.G. Sayegh, J.H. Vera, Lattice - Model Expressions for the Combinatorial Entropy of Liquid Mixtures: A Critical Discussion, The Chemical Engmeerrng Journal,. 19 (1980) 1 10. [114] G.M. Kontogeorgis, G.K. Folas, Chapter 9, Appendix A, in: Thermodynamic Models for Industrial Applications: From Classical and Advanced Mixing Rules to Association Theories, John Wiley & Sons, 2009. h ttps://www.wiley.com/go/kontogeorgis. [115] J. Gross, G. Sadowski, Application of the Perturbed - Chain SAFT Equation of State to Associating Systems, Ind. Eng. Chem. Res. 41 (2002) 5510 5515. doi:10.1021/ie010954d. [116] N.F. Carnahan, K.E. Starling, Equati on of State for Nonattracting Rigid Spheres, The Journal of Chemical Physics. 51 (1969) 635 636. doi:10.1063/1.1672048. [117] L.S. Budantseva, T.M. Lesteva, M.S. Nemtsov, Liquid - vapor equilibrium in methanol - C6 hydrocarbons of different classes, Zh. Fiz. K him. 49 (1975) 260. [118] H. Katayama, M. Ichikawa, Liquid - liquid equilibria of three ternary systems: methanol - heptane including 1, 3 - dioxolane, 1, 4 - dioxane and tetrahydropyran in the range of 253.15 to 303.15 K, Journal of Chemical Engineering of Japan. 28 (1995) 412 418. 203 [119] J. Gmehling, B. Krentscher, Excess enthalpies of 12 binary liquid mixtures containing cyclohexane at elevated temperatures and pressures (up to 416 K and 1.9 MPa), ELDATA Int. Electron. J. Phys. - Chem. Data. 1 (1995) 181 190. [120] R.A. Wilsak, S.W. Campbell, G. Thodos, Vapor liquid equilibrium measurements for the n - pentane methanol system at 372.7, 397.7 and 422.6 K, Fluid Phase Equilibria. 33 (1987) 157 171. [121] Fu, J., Wang, K., Hu, Y., Studies on vapor - liquid and liquid - liqui d vapor equilibria for the ternary system methanol - methyl methacrylate - water (I) three binary systems, Chinese Journal of Chemical Engineering. 4 (1989) 1 13. [122] C. Yang, X. Yin, S. Ma, Organic salt effect of tetramethylammonium bicarbonate on the vapor liquid equilibrium of the dimethyl carbonate+ methanol system, Journal of Chemical & Engineering Data. 57 (2011) 66 71. [ 123 ] A. Arce, J. Mart nez - Ageitos, E. Rodil, A. Soto, Measurement and prediction of isobaric vapour liquid equilibrium data of the sy stem ethanol+ methanol+ 2 - methoxy - 2 - methylpropane, Fluid Phase Equilibria. 146 (1998) 139 153. [124] Z.S. Kooner, D.V. Fenby, Vapour pressure study of the deuterium exchange reaction in methanol - ethanol systems: equilibrium constant determination, Australi an Journal of Chemistry. 33 (1980) 1943 1946. [125] K.L. Butcher, W.I. Robinson, An apparatus for determining high - pressure liquid - vapour equilibrium data. I. The methanol ethanol system, Journal of Chemical Technology and Biotechnology. 16 (1966) 289 292. [126] H.D. Pflug, A.E. Pope, G.C. Benson, Heats of mixing of normal alcohols at 25. deg., Journal of Chemical & Engineering Data. 13 (1968) 408 410. [127] H.S. Myers, VL equilibrium for naphthenes and paraffins, Petroleum Refiner. 36 (1957) 175 178. [128] J.M. Sturtevant, P.A. Lyons, Application of flow calorimetry to the determination of the enthalpies of mixing of organic liquids, The Journal of Chemical Thermodynamics. 1 (1969) 201 209. [129] C. Narasigadu, K. Moodley, J.D. Raal, Vapor Liquid Equilibriu m Measurements of Ethanol and n - Nonane or n - Decane Binary Mixtures with Large Relative Volatility, Journal of Chemical & Engineering Data. 63 (2018) 1240 1248. [130] L. Wang, G.C. Benson, B.C. - Y. Lu, Excess enthalpies of (ethanol+ hexane+ decane or dodecan e) at the temperature 298.15 K, The Journal of Chemical Thermodynamics. 24 (1992) 1135 1143. [131] K.S. Yuan, J.C.K. Ho, A.K. Koshpande, B.C. - Y. Lu, Vapor - Liquid Equilibria., J. Chem. Eng. Data. 8 (1963) 549 559. doi:10.1021/je60019a024. 204 [132] I. Nagata, K . Kazuma, Heats of mixing for the ternary system ethanol - 1 - propanol - cyclohexane at 25. degree. C, Journal of Chemical and Engineering Data. 22 (1977) 79 84. [133] E. Tusel - Langer, J.M.G. Alonso, M.A.V. Olfos, R.N. Lichtenthaler, Excess enthalpies of mixtur es containing n - heptane, methanol and methyl tert - butyl ether (MTBE), J Solution Chem. 20 (1991) 153 163. doi:10.1007/BF00651647. [134] P.M. Mathias, H.Z. Kister, Effect of Phase - Equilibrium Uncertainties on Ethyl Acetate Purification, J. Chem. Eng. Data. 62 (2017) 2872 2883. doi:10.1021/acs.jced.7b00172. [135] V. Flemr, Note on excess Gibbs energy equations based on local composition concepts., Collection of Czechoslovak Chemical Communications. 41 (1976) 3347 3349. [136] S. - T. Lin, M. - K. Hsieh, C. - M. Hsie h, C. - C. Hsu, Towards the development of theoretically correct liquid activity coefficient models, The Journal of Chemical Thermodynamics. 41 (2009) 1145 1153. doi:10.1016/j.jct.2009.05.002. [137] N. Asprion, H. Hasse, G. Maurer, FT - IR spectroscopic invest igations of hydrogen bonding in alcohol - hydrocarbon solutions, Fluid Phase Equilibria. 186 (2001) 1 25. doi:10.1016/s0378 - 3812(01)00363 - 6. [138] N. Asprion, H. Hasse, G. Maurer, Thermodynamic and IR spectroscopic studies of solutions with simultaneous asso ciation and solvation, Fluid Phase Equilibria. 208 (2003) 23 51. doi:10.1016/S0378 - 3812(02)00317 - 5. [139] S.D. Williams, T.J. Johnson, S.W. Sharpe, V. Yavelak, R.P. Oates, C.S. Brauer, Quantitative vapor - phase IR intensities and DFT computations to predict absolute IR spectra based on molecular structure: I. Alkanes, Journal of Quantitative Spectroscopy and Radiative Transfer. 129 (2013) 298 307. doi:10.1016/j.jqsrt.2013.07.005. [140] N. Kumar, A. Bansal, G.S. Sarma, R.K. Rawal, Chemometrics tools used in a nalytical chemistry: an overview, Talanta. 123 (2014) 186 199. doi:10.1016/j.talanta.2014.02.003. [141] R.N. Jones, The intensities of the infra - red absorption bands of n - paraffin hydrocarbons, Spectrochimica Acta. 9 (1957) 235 251. doi:10.1016/0371 - 1951(5 7)80137 - 4. [142] A.S. Wexler, Infrared determination of structural units in organic compounds by integrated intensity measurements: Alkanes, alkenes and monosubstituted alkyl benzenes, Spectrochimica Acta. 21 (1965) 1725 1742. doi:10.1016/0371 - 1951(65)8008 5 - 6. [143] by infrared spectroscopy. II. Variation of the intensity of the antisymmetrical stretching vibration bands of C - H bonds in CH 3 and CH 2 groups, Collect. Czech. Che m. Commun., CCCC. 32 (1967) 1903 1912. doi:10.1135/cccc19671903. [144] K.M. Murdoch, T.D. Ferris, J.C. Wright, T.C. Farrar, Infrared spectroscopy of ethanol clusters in ethanol hexane binary solutions, The Journal of Chemical Physics. 116 (2002) 5717 5724. doi:10.1063/1.1458931. 205 [145] E.E. Tucker, E.D. Becker, Alcohol association studies. II. Vapor pressure, 220 - MHz proton magnetic resonance, and infrared investigations of tert - butyl alcohol association in hexadecane, The Journal of Physical Chemistry. 77 ( 2012) 1783 1795. doi:10.1021/j100633a012. [146] X. Wu, Y. Chen, T. Yamaguchi, Hydrogen bonding in methanol studied by infrared spectroscopy, Journal of Molecular Spectroscopy. 246 (2007) 187 191. doi:10.1016/j.jms.2007.09.012. [147] J.T. Reilly, A. Thomas, A.R. Gibson, C.Y. Luebehusen, M.D. Donohue, Analysis of the Self - Association of Aliphatic Alcohols Using Fourier Transform Infrared (FT - IR) Spectroscopy, Ind. Eng. Chem. Res. 52 (2013) 14456 14462. doi:10.1021/ie302174r. [148] K. Ohta, K. Tominaga, Vibrat ional population relaxation of hydrogen - bonded phenol complexes in solution: Investigation by ultrafast infrared pump probe spectroscopy, Chemical Physics. 341 (2007) 310 319. doi:10.1016/j.chemphys.2007.07.025. [149] K. Mazur, M. Bonn, J. Hunger, Hydrogen Bond Dynamics in Primary Alcohols: A Femtosecond Infrared Study, J. Phys. Chem. B. 119 (2015) 1558 1566. doi:10.1021/jp509816q. [150] A. Choperena, P. Painter, An infrared spectroscopic study of hydrogen bonding in ethyl phenol: A model system for polymer phenolics, Vibrational Spectroscopy. 51 (2009) 110 118. doi:10.1016/j.vibspec.2008.11.008. [151] A. Choperena, P. Painter, Hydrogen Bonding in Polymers: Effect of Temperature on the OH Stretching Bands of Poly(vinylphenol), Macromolecules. 42 (2009) 6159 6165. doi:10.1021/ma900928z. [152 ] N. Sheppard, Infrared Spectroscopy and Hydrogen Bonding - Band - Widths and Frequency Hydrogen Bonding Held at Ljubljana, 29 July 3 August 1957, Pergamon, 1959: pp. 85 105. doi:10.1016/B978 - 0 - 08 - 009140 - 2.50013 - 0. [153] J. - M. Andanson, J. - C. Soetens, T. Tassaing, M. Besnard, Hydrogen bonding in supercritical tert - butanol assessed by vibrational spectroscopies and molecular - dynamics simulations., The Journal of Chemical Physics . 122 (2005) 174512. doi:10.1063/1.1886730. [154] I. Doroshenko, V. Pogorelov, V. Sablinskas, V. Balevicius, Matrix - isolation study of cluster formation in methanol: O H stretching region, Journal of Molecular Liquids. 157 (2010) 142 145. doi:10.1016/j.mol liq.2010.09.003. [155] D. Wandschneider, M. Michalik, A. Heintz, Spectroscopic and thermodynamic studies of liquid n - butanol+n - hexane and +cyclohexane mixtures based on quantum mechanical ab initio calculations of n - butanol clusters, Journal of Molecular L iquids. 125 (2006) 2 13. doi:10.1016/j.molliq.2005.11.011. 206 [156] P.E. Hansen, J. Spanget - Larsen, On prediction of OH stretching frequencies in intramolecularly hydrogen bonded systems, Journal of Molecular Structure. 1018 (2012) 8 13. doi:10.1016/j.molstru c.2012.01.011. [157] M. Koné, B. Illien, C. Laurence, J. Graton, Can Quantum - Mechanical Calculations Yield Reasonable Estimates of Hydrogen - Bonding Acceptor Strength? The Case of Hydrogen - Bonded Complexes of Methanol, The Journal of Physical Chemistry A. 1 15 (2011) 13975 13985. doi:10.1021/jp209200w. [158] K. Ohno, T. Shimoaka, N. Akai, Y. Katsumoto, Relationship between the broad OH stretching band of methanol and hydrogen - bonding patterns in the liquid phase., The Journal of Physical Chemistry. A. 112 (20 08) 7342 8. doi:10.1021/jp800995m. [159] J. Spanget - Larsen, B.K.V. Hansen, P.E. Hansen, OH stretching frequencies in systems with intramolecular hydrogen bonds: Harmonic and anharmonic analyses, Chemical Physics. 389 (2011) 107 115. doi:10.1016/j.chemphys. 2011.09.011. [160] A. Staib, D. Borgis, A quantum multi - mode molecular dynamics approach to the vibrational spectroscopy of solvated hydrogen - bonded complexes, Chemical Physics Letters. 271 (1997) 232 240. doi:10.1016/S0009 - 2614(97)00470 - 3. [161] M.K. Ghos h, J. Lee, C.H. Choi, M. Cho, Direct Simulations of Anharmonic Infrared Spectra Using Quantum Mechanical/Effective Fragment Potential Molecular Dynamics (QM/EFP - MD): Methanol in Water, J. Phys. Chem. A. 116 (2012) 8965 8971. doi:10.1021/jp306807v. [162] J. Wang, R.J. Boyd, A. Laaksonen, A hybrid quantum mechanical force field molecular dynamics simulation of liquid methanol: Vibrational frequency shifts as a probe of the quantum mechanical/molecular mechanical coupling, The Journal of Chemical Physics. 104 (1996) 7261 7269. doi:10.1063/1.471439. [163] S.A. Corcelli, J.L. Skinner, Infrared and Raman Line Shapes of Dilute HOD in Liquid H2O and D2O from 10 to 90 °C, J. Phys. Chem. A. 109 (2005) 6154 6165. doi:10.1021/jp0506540. [164] S.A. Corcelli, C.P. Lawrenc e, J.L. Skinner, Combined electronic structure/molecular dynamics approach for ultrafast infrared spectroscopy of dilute HOD in liquid H2O and D2O, The Journal of Chemical Physics. 120 (2004) 8107 8117. doi:10.1063/1.1683072. [165] J.R. Schmidt, S.A. Corce lli, J.L. Skinner, Pronounced non - Condon effects in the ultrafast infrared spectroscopy of water, The Journal of Chemical Physics. 123 (2005) 044513. doi:10.1063/1.1961472. [166] B. Auer, R. Kumar, J.R. Schmidt, J.L. Skinner, Hydrogen Bonding and Raman, IR , and 2D - Sciences of the United States of America. 104 (2007) 14215 14220. 207 [167] S.M. Gruenbaum, C.J. Tainter, L. Shi, Y. Ni, J.L. Skinner, Robustness of Frequency, Transit ion Dipole, and Coupling Maps for Water Vibrational Spectroscopy, J. Chem. Theory Comput. 9 (2013) 3109 3117. doi:10.1021/ct400292q. [168] Mode in Alcohols, J. Phys. Chem. A. 121 (2017) 5823 5833. doi:10.1021/acs.jpca.7b05836. [169] K.B. Whetsel, J.. Lady, Self - association of phenol in nonpolar solvents, in: R.A. Friedel (Ed.), Spectrometry of Fuels, Plenum Press, New York, 1970. [170] A. Hall, J.L. Wood, Origin of the structu re of the OH stretching band of some phenols in solution, Spectrochimica Acta Part A: Molecular Spectroscopy. 23 (1967) 2657 2667. doi:10.1016/0584 - 8539(67)80157 - 0. [171] S. Woutersen, U. Emmerichs, H.J. Bakker, A femtosecond midinfrared pump probe study o f hydrogen - bonding in ethanol, The Journal of Chemical Physics. 107 (1997) 1483 1490. doi:10.1063/1.474501. [172] K. Shinokita, A.V. Cunha, T.L.C. Jansen, M.S. Pshenichnikov, Hydrogen bond dynamics in bulk alcohols, The Journal of Chemical Physics. 142 (20 15) 212450. doi:10.1063/1.4921574. [173] Journal of Chemical Physics. 22 (1954) 1697 1701. doi:10.1063/1.1739878. [174] D. Deng, H. Li, J. Yao, S. Han, Simple local composition model for 1H NMR chemical shift of mixtures, Chemical Physics Letters. 376 (2003) 125 129. doi:10.1016/S0009 - 2614(03)00996 - 5. [175] D.A. Case, V. Babin, J.T. Berryman, R.M. Betz, Q. Cai, D.S. Cerutti, T.E. Cheatham, III, T.A. Darden, R.E., Duke, H. Gohlke, A.W. Goetz, S. Gusarov, N. Homeyer, P. Janowski, J. Kaus, I. Kolossváry, A. Kovalenko, T.S. Lee, S. LeGrand, T. Luchko, R. Luo, B. Madej, K.M. Merz, F. Paesani, D.R. Roe, A. Roitberg, C. Sagui, R. Salomon - Ferrer, G. Seabra, C.L. Simmerling, W. Smith, J. S wails, R.C. Walker, J. Wang, R.M. Wolf, X., Wu and P.A. Kollman (2014), AMBER 14, University of California, San Francisco, n.d. [176] L. Martínez, R. Andrade, E.G. Birgin, J.M. Martínez, PACKMOL: A package for building initial configurations for molecular dynamics simulations, Journal of Computational Chemistry. 30 (n.d.) 2157 2164. doi:10.1002/jcc.21224. [177] M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, G. A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H. P. Hratchian, A. F. Izmaylov, J. B loino, G. Zheng, J. L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. Bearpark, J. J. Heyd, E. Brothers, K. N. Kudin, V . N. Staroverov, T. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J. M. Millam, M. Klene, J. E. Knox, J. B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R. E. Stratmann, O . Yazyev, A. J. Austin, R. Cammi, C. Pomelli, J. W. Ochterski, R. L. 208 Martin, K. Morokuma, V. G. Zakrzewski, G. A. Voth, P. Salvador, J. J. Dannenberg, S. Dapprich, A. D. Daniels, O. Farkas, J. B. Foresman, J. V. Ortiz, J. Cioslowski,, and D. J. Fox, Gaussi an 09, Revision D.01, Gaussian, Inc., Wallingford CT, 2013. [178] A.D. Becke, Density - functional thermochemistry. V. Systematic optimization of exchange - correlation functionals, The Journal of Chemical Physics. 107 (1997) 8554 8560. doi:10.1063/1.475007. [ 179] R. Zheng, Y. Sun, Q. Shi, Theoretical study of the infrared and Raman line shapes of liquid methanol, Phys. Chem. Chem. Phys. 13 (2011) 2027 2035. doi:10.1039/C0CP01145B. [180] H. Torii, Time - Domain Calculations of the Infrared and Polarized Raman Spe ctra of Tetraalanine in Aqueous Solution, J. Phys. Chem. B. 111 (2007) 5434 5444. doi:10.1021/jp070301w. [181] Killian W.G.; Unpublished experiments and calculations from Lira Lab, Michigan State University; Personal Communication, 2018, (n.d.). [182] J. S panget - Output, (2015). doi:10.13140/RG.2.1.4181.6160. [183] NIST Computational Chemistry Comparison and Benchmark Database; NIST Standard Reference Database Number 101; Release 19, A pril 2018, Editor: Russell D. Johnson III; http://cccbdb.nist.gov/; DOI:10.18434/T47C7Z, n.d. [184] C. Van Ness, B. Hollingerlb, Infrared Spectra and the Thermodynamics of Alcohol - Hydrocarbon Systems, The Journal of Physical Chemistry. 71 (1967) 1483 1494. [185] A.N. Fletcher, C.A. Heller, Self - association of alcohols in nonpolar solvents, J. Phys. Chem. 71 (1967) 3742 3756. doi:10.1021/j100871a005. [186] K. Shinomiya, T. Shinomiya, An Equilibrium Model of the Self - Association of 1 - and 3 - Pentanols in Hept ane, The Chemical Society of Japan. 63 (1990) 1093 1097. [187] C.M. Huggins, G.C. Pimentel, Systematics of the Infrared Spectral Properties of Hydrogen Bonding Systems: Frequency Shift, Half Width and Intensity, J. Phys. Chem. 60 (1956) 1615 1619. doi:10.1 021/j150546a004. [188] T. Ebukuro, A. Takami, Y. Oshima, S. Koda, Raman spectroscopic studies on hydrogen bonding in methanol and methanol/water mixtures under high temperature and pressure, The Journal of Supercritical Fluids. 15 (1999) 73 78. doi:10.1016 /S0896 - 8446(98)00126 - 0. [189] - Larsen, Intramolecular hydrogen bonding in myricetin and myricitrin. Quantum chemical calculations and vibrational spectroscopy, Journal of Molecular Structure. 1131 (2017) 242 249. doi:10.1016/j.molst ruc.2016.11.069. 209 [190] F. Palombo, M. Paolantoni, P. Sassi, A. Morresi, R.S. Cataliotti, Spectroscopic studies of 139 146. [191] T. Yamaguchi, C.J. Benmore, A.K. So per, The structure of subcritical and supercritical methanol by neutron diffraction, empirical potential structure refinement, and spherical harmonic analysis, The Journal of Chemical Physics. 112 (2000) 8976 8987. doi:10.1063/1.481530. [192] T. Yamaguchi, K. Hidaka, A.K. Soper, The structure of liquid methanol revisited: a neutron 1168. doi:10.1080/00268979909483060. [193] C.J. Benmore, Y.L. Loh, The structure of liquid ethanol: A neutron diffraction and molecular dynamics study, The Journal of Chemical Physics. 112 (2000) 5877 5883. doi:10.1063/1.481160. [194] C. Czeslik, J. Jonas, Pressure and temperature dependence of hydrogen - bond strength in methanol clusters, Chemical Physi cs Letters. 302 (1999) 633 638. doi:10.1016/S0009 - 2614(99)00170 - 0. [195] G.C. Pimentel, A.L. McClellan, The hydrogen bond, 2nd Edition, W.H. Freeman, San Francisco, 1960. [196] P. Schuster, G theory and experiments, North - [197] S.N. Vinogradov, R.H. Linnell, Hydrogen bonding, Van Nostrand Reinhold, New York, 1971. [198 ] D.S. Bulgarevich, K. Otake, T. Sako, T. Sugeta, Y. Takebayashi, C. Kamizawa, D. Shintani, A. Negishi, C. Tsurumi, Hydrogen bonding in supercritical methanol studied by infrared spectroscopy, The Journal of Chemical Physics. 116 (2002) 1995 2003. doi:10.1 063/1.1431585. [199] R. Böhmer, C. Gainaru, R. Richert, Structure and dynamics of monohydroxy alcohols Milestones towards their microscopic understanding, 100 years after Debye, Physics Reports. 545 (2014) 125 195. doi:10.1016/j.physrep.2014.07.005. [200] S.J. Barlow, G.V. Bondarenko, Y.E. Gorbaty, T. Yamaguchi, M. Poliakoff, An IR Study of Hydrogen Bonding in Liquid and Supercritical Alcohols, J. Phys. Chem. A. 106 (2002) 10452 10460. doi:10.1021/jp0135095. [201] O. Gereben, L. Pusztai, Investigation of th e Structure of Ethanol Water Mixtures by Molecular Dynamics Simulation I: Analyses Concerning the Hydrogen - Bonded Pairs, J. Phys. Chem. B. 119 (2015) 3070 3084. doi:10.1021/jp510490y. 210 [202] ure of Ethanol Water Mixtures by Molecular Dynamics Simulation I: Analyses Concerning the Hydrogen - Bonded 4564. doi:10.1021/acs.jpcb.5b02167. [203] O. Gereben, Ring structure analysis of ethanol water mixtures, Jou rnal of Molecular Liquids. 211 (2015) 812 820. doi:10.1016/j.molliq.2015.08.010. [204] Y. Vaskivskyi, I. Doroshenko, Y. Chernolevska, V. Pogorelov, G. Pitsevich, Spectroscopic studies of clusterization of methanol molecules isolated in a nitrogen matri x, Low Temperature Physics. 43 (2017) 1415 1419. doi:10.1063/1.5012794. [205] R. Laenen, C. Rauscher, A structural model for associated liquid ethanol developed from transient spectroscopy, The Journal of Chemical Physics. 107 (1997) 9759 9763. doi:10.1063 /1.475273. [206] T. Häber, U. Schmitt, M.A. Suhm, FTIR - spectroscopy of molecular clusters in pulsed supersonic slit - jet expansions, Physical Chemistry Chemical Physics. 1 (1999) 5573 5582. doi:10.1039/A907264K. [207] M. Paolantoni, P. Sassi, A. Morresi, R. S. Cataliotti, Infrared study of 1 - octanol liquid structure, Chemical Physics. 310 (2005) 169 178. doi:10.1016/j.chemphys.2004.10.027. [208] P. Gans, Data fitting in the chemical sciences: by the method of least squares, Wiley, Chichester, 1992. [209] J. - S . Chen, K. - T. Yeh, D. - Y. Kao, J.K. Baird, Toward a Better Determination of Infrared Molar Absorptivities and Dimer Formation Constants in Self - association Through Hydrogen Bonding: 3 - Ethyl - 2 - methyl - 3 - pentanol in Tetrachloroethylene as an Example, Journal o f the Chinese Chemical Society. 64 (n.d.) 1156 1163. doi:10.1002/jccs.201600865. [210] H. Kempter, R. Mecke, Spektroskopische Bestimmung von Assoziations - Gleichgewichten, Die Naturwissenschaften. 27 (1939) 583 584. doi:10.1007/BF01496162. [211] G.D. Ikonom ou, M.D. Donohue, Thermodynamics of hydrogen - bonded molecules: The associated perturbed anisotropic chain theory, AIChE Journal. 32 (n.d.) 1716 1725. doi:10.1002/aic.690321015. [212] G.D. Ikonomou, M.D. Donohue, Compact: a simple equation of state for asso ciated molecules, Fluid Phase Equilibria. 33 (1987) 61 90. doi:10.1016/0378 - 3812(87)87004 - 8. [213] J.M. Bakke, I. Hampson, S. Kumar, J. Gallagher, R. Datema, J. Chattopadhyaya, Hydrogen Bonding of t - Butyl Alcohol at Low Temperatures and Concentrations., Ac ta Chemica Scandinavica. 39b (1985) 15 23. doi:10.3891/acta.chem.scand.39b - 0015. [214] G. Brink, L. Glasser, Dielectric studies of molecular association. A model for the association of ethanol in dilute solution, J. Phys. Chem. 82 (1978) 1000 1005. doi:10. 1021/j100498a007. 211 [215] C.A. Eckert, M.M. McNiel, B.A. Scott, L.A. Halas, NMR Measurements of chemical theory equilibrium constants for hydrogen - bonded solutions, AIChE Journal. 32 (1986) 820 828. doi:10.1002/aic.690320512. [216] H.S. Gutowsky, A. Saika, D issociation, Chemical Exchange, and the Proton Magnetic Resonance in Some Aqueous Electrolytes, The Journal of Chemical Physics. 21 (1953) 1688 1694. doi:10.1063/1.1698644.