Modelado-termodinamico-y-simulacion-molecular-de-mezclas-ternarias-agua-hidrocarburo-tensoactivo

Vista previa del material en texto

UNIVERSIDAD NACIONAL AUTÓNOMA DE MEXICO 
 
 
PROGRAMA DE MAESTRÍA Y DOCTORADO EN CIENCIAS QUÍMICAS 
 
 
 
 
MODELADO TERMODINÁMICO Y SIMULACIÓN MOLECULAR DE MEZCLAS 
TERNARIAS AGUA–HIDROCARBURO–TENSOACTIVO 
 
 
 
 
TESIS 
 
PARA OPTAR POR EL GRADO DE 
 
 
DOCTOR EN CIENCIAS 
 
 
 
 
PRESENTA 
 
M. en I. EDGAR GALICIA ANDRÉS 
 
 
 
DIRECTOR DE TESIS: 
DR. MILTON THADEU GARCÍA MEDEIROS DE OLIVEIRA 
FACULTAD DE QUÍMICA 
 
 
 
CIUDAD UNIVERSITARIA, CD. MX., DICIEMBRE 2017 
 
 
UNAM – Dirección General de Bibliotecas 
Tesis Digitales 
Restricciones de uso 
 
DERECHOS RESERVADOS © 
PROHIBIDA SU REPRODUCCIÓN TOTAL O PARCIAL 
 
Todo el material contenido en esta tesis esta protegido por la Ley Federal 
del Derecho de Autor (LFDA) de los Estados Unidos Mexicanos (México). 
El uso de imágenes, fragmentos de videos, y demás material que sea 
objeto de protección de los derechos de autor, será exclusivamente para 
fines educativos e informativos y deberá citar la fuente donde la obtuvo 
mencionando el autor o autores. Cualquier uso distinto como el lucro, 
reproducción, edición o modificación, será perseguido y sancionado por el 
respectivo titular de los Derechos de Autor. 
 
 
 
2 
 
 
 
UNIVERSIDAD NACIONAL AUTÓNOMA DE MÉXICO 
 
 
 
PROGRAMA DE MAESTRÍA Y DOCTORADO EN CIENCIAS 
QUÍMICAS 
 
 
 
 
MODELADO TERMODINÁMICO Y SIMULACIÓN MOLECULAR DE 
MEZCLAS TERNARIAS AGUA–HIDROCARBURO–TENSOACTIVO 
 
 
 
T E S I S 
PARA OPTAR POR EL GRADO DE 
 
 
DOCTOR EN CIENCIAS 
 
 
P R E S E N T A 
M. en I. EDGAR GALICIA ANDRÉS 
 
DIRECTOR DE TESIS 
DR. MILTON THADEU GARCÍA MEDEIROS DE OLIVEIRA 
FACULTAD DE QUÍMICA 
 
 Ciudad de México, Diciembre 2017 
 
3 
 
 
Contenido 
Contenido................................................................................................................................ 3 
Resumen ................................................................................................................................. 8 
Abstract ................................................................................................................................... 9 
1 Introducción .................................................................................................................. 10 
1.1 Plan de tesis............................................................................................................ 13 
Objetivos ............................................................................................................................... 16 
Hipótesis ............................................................................................................................... 17 
2 Marco teórico ................................................................................................................ 18 
2.1 Dinámica molecular ............................................................................................... 19 
2.1.1 Ecuaciones de movimiento ............................................................................. 20 
2.1.2 Algoritmo salto de rana .................................................................................. 20 
2.1.3 Condiciones de frontera .................................................................................. 21 
2.1.4 Interacciones Coulombicas (Sumas de Ewald) .............................................. 23 
2.1.5 Funciones potenciales ..................................................................................... 24 
2.2 Propiedades interfaciales ....................................................................................... 30 
2.2.1 Perfiles de densidad ........................................................................................ 31 
2.2.2 Tensión superficial .......................................................................................... 31 
2.2.3 Parámetro de orden orientacional ................................................................... 32 
2.2.4 Distancia de extremo a extremo y radio de giro ............................................. 33 
2.2.5 Coeficiente de difusión ................................................................................... 34 
2.2.6 Energía libre.................................................................................................... 35 
2.3 Teoría gradiente ..................................................................................................... 35 
2.4 Ecuaciones de estado ............................................................................................. 36 
2.4.1 Teoría estadística de fluidos asociantes (SAFT) ............................................. 37 
2.4.2 Ecuación de estado SRK................................................................................. 40 
3 Metodología .................................................................................................................. 44 
3.1 Detalles de la simulación ....................................................................................... 44 
4 Resultados y discusión.................................................................................................. 48 
4.1 n–Alcanos .............................................................................................................. 48 
4 
 
4.1.1 Perfiles de densidad ........................................................................................ 48 
4.1.2 Tensión superficial .......................................................................................... 51 
4.1.3 Parámetro de orden orientacional ................................................................... 53 
4.1.4 Distancia de extremo a extremo ..................................................................... 55 
4.1.5 Radio de giro .................................................................................................. 56 
4.1.6 Coeficiente de difusión ................................................................................... 57 
4.1.7 Perfiles de energía libre .................................................................................. 59 
4.2 Agua–Hidrocarburos .............................................................................................. 60 
4.2.1 Densidades y tensión interfacial ..................................................................... 61 
4.2.2 Parámetros de orden orientacional ................................................................. 63 
4.2.3 Distancia de extremo a extremo y radios de giro ........................................... 68 
4.2.4 Coeficientes de difusión ................................................................................. 69 
4.2.5 Tensoactivos ................................................................................................... 70 
4.3 Agua–hidrocarburos–tensoactivos ......................................................................... 74 
4.3.1 Tensión interfacial y perfil de densidad .......................................................... 75 
4.3.2 Parámetro de orden orientacional ................................................................... 78 
4.3.3 Coeficientes de difusión ................................................................................. 81 
5 Conclusiones ................................................................................................................. 84 
Trabajos a futuro ................................................................................................................... 86 
 
 
5 
 
Jurado asignado: 
Presidente: Dra. Jacqueline Quintana 
Vocal: Dr. José Alejandre Ramírez 
Vocal: Dr. Jesús Gracia Fadrique 
Vocal: Dr. Jorge Balmaseda Era 
Secretario: Dr. David Philip Sanders 
 
Lugar donde se realizó el trabajo: 
Laboratorio 311, Departamento de Fisicoquímica, Facultad de Química, UNAM. 
 
Nombre del Tutor: 
Dr. Milton Thadeu García Medeiros de Oliveira 
 
6 
 
AgradecimientosA la Universidad Nacional Autónoma de México por brindar los recursos para llevar a cabo 
esta investigación. 
Al Dr. Milton Medeiros por permitirme formar parte de su proyecto y brindarme su apoyo y 
amistad. 
Al CONACYT por otorgar la beca del periodo 2014-2 – 2017-1 a nombre de Edgar Galicia 
Andrés (No. De Becario: 255431). 
 
A los miembros del jurado por dedicar su tiempo y sus valiosos comentarios: 
Dra. Jacqueline Quintana Hinojosa 
Dr. José Alejandre Ramírez 
Dr. Jesús Gracia Fadrique 
Dr. Jorge Balmaseda Era 
Dr. David Philip Sanders 
 
Al Dr. Joaquín Barroso Flores, miembro del comité tutor, por sus valiosos comentarios para 
enriquecer el proyecto. 
 
Publicación: 
E. Galicia–Andrés, M. Medeiros, “Vapour–Liquid Interfacial Properties of n–alkanes” 
Journal of Molecular Liquids, Vol. 248 (2017) 253–263. 
 
7 
 
Dedicatoria 
A Dios por todas las bendiciones que me ha dado. 
A mi esposa MoNsE por todo su amor, cariño, paciencia por ser una fuente de inspiración y 
sobre todo, su apoyo. ¡Gracias princesa! 
A mis padres por todo el amor y apoyo que me han dado. 
 
8 
 
Resumen 
El proyecto consiste en el cálculo numérico de propiedades interfaciales de mezclas de agua 
con n–alcanos y tensoactivos de la familia de alcoholes etoxilados usando simulaciones 
moleculares, así como elementos de la teoría gradiente de interfases fluidas. 
Los sistemas de interés comprenden el equilibrio líquido–vapor, ELV, de sustancias puras y 
tensoactivos en solución acuosa, el equilibrio líquido–líquido, ELL, de mezclas agua–
hidrocarburo y agua–hidrocarburo–tensoactivos. Los alcanos lineales estudiados son n–
hexano, n–heptano, n–octano, n–nonano y n–decano, mientras que los tensoactivos son 
C12E2, C12E3, C12E4 y C12E5. 
De los resultados obtenidos se destaca la orientación preferencial de los n–alcanos, del agua 
y de los tensoactivos en diferentes regiones de la interfase, así como el cambio de las 
propiedades dinámicas a lo largo de la interfase. 
Los n–alcanos tienden a orientarse perpendiculares a la interfase, cercanos al bulto de la fase 
líquido y vapor, mientras que, en la superficie divisoria de Gibbs (SDG) se orientan paralelos 
a la interfase. Esta orientación se ve afectada por la presencia de moléculas de agua, las cuales 
incrementan la intensidad de orientación de los alcanos en la SDG y disminuye el número de 
moléculas con orientación perpendicular en la cercanía de la fase orgánica y acuosa. 
Por otro lado, las moléculas de agua tienden a orientar su vector OH para exponer su átomo 
de hidrogeno al bulto de la fase líquida y vapor. Similar a los alcanos, las moléculas de agua 
orientan su vector de momento dipolar paralelo a la interfase en la SDG. El agua no se ve 
afectada por la presencia de hidrocarburos; sin embargo, la presencia de tensoactivos cambia 
la orientación preferencial tanto del agua como de los n–alcanos, haciendo que los alcanos 
se orienten perpendiculares a la interfase en equilibrio con una fase acuosa y el agua cambia 
la orientación del vector OH en la cercanía de la fase orgánica. 
Con respecto a la difusividad lateral, se observa una alta movilidad en la región interfacial 
que se ve alterada por la presencia de tensoactivo, siendo un impedimento para la movilidad 
entre dos fases. 
 
 
9 
 
Abstract 
The project consists of the calculation of interfacial properties of water, n–alkanes and 
surfactants from the alcohol ethoxylates family mixtures using molecular simulations and 
elements from square gradient theory. 
The systems of interest include the vapor–liquid equilibrium, VLE, of pure substances and 
surfactants in solution, liquid–liquid equilibrium, LLE, of water–alkane mixtures and water–
alkane–surfactants. The n–alkanes under study are n–hexane, n–heptane, n–octane, n–
nonane, n–decane, whilst the surfactants are C12E2, C12E3, C12E4 y C12E5. 
From the results obtained, we highlight the orientational preference of n–alkanes, water and 
surfactants in different regions of the interface, as well as the change of the dynamic 
properties along the interface. 
The n–alkanes tend to orient perpendicular to the interface, close to the liquid and vapour 
bulk phases, while at the Gibbs dividing surface (GDS) the n–alkanes orient parallel to the 
interface. The orientation is affected by the presence of water molecules, which increase the 
orientation intensity of the alkanes at the GDS and decrease the number of molecules with 
perpendicular orientation close to the organic and aqueous phases. 
On the other hand, water molecules tend to orient their OH vector to expose their hydrogen 
atoms to the liquid and vapour bulk phase. Similar to the alkanes, water molecules orient 
their dipole moment vector parallel to the interface at the SDG. Water is not affected by the 
presence of hydrocarbons; meanwhile, the presence of surfactants changes the preferential 
orientations of water and n–alkane molecules, making the alkanes orient themselves 
perpendicular to the interface in equilibrium with water, while water molecules change the 
orientation of OH vector close to the organic phase. 
Regarding lateral diffusion, a high mobility is observed in the interfacial region, which is 
altered by the presence of the surfactants, being an impediment for the mobility between the 
two phases. 
 
10 
 
1 Introducción 
Las propiedades interfaciales son de gran importancia en diversos campos de la industria 
(energéticos, alimentos, cosméticos, etc.) para la separación de fases, emulsificación, 
extracción de sustancias de interés o incluso la administración de fármacos, por citar algunos 
ejemplos. Dentro de los procesos para la separación de fases, la emulsificación es uno de los 
más socorridos, con el uso de los denominados tensoactivos que tienen la capacidad de abatir 
la tensión interfacial. 
En la industria del petroleo, existen diferentes etapas que comprenden el proceso de 
extracción en los yacimientos de petróleo [1]. La recuperación primaria, que consiste en 
extraer los hidrocarburos con la misma presión del yacimiento; La recuperación secundaria 
consiste en la inyección de un fluido (agua, reinyección de gas natural u otros gases al fondo) 
para incrementar la presión del pozo y reducir la densidad promedio en el pozo. Este proceso 
recupera alrededor del 30% [2]; La recuperación asistida de pertróleo, o recuperación 
terciaria son métodos que reduciden la viscosidad, haciendo más fácil su extracción. Es 
común la inyección de polímeros para inundar posteriormente con agua, surfactantes para 
reducir la tensión interfacial y movilizar el petróleo a los puntos de extracción [1,3]. Debido 
a esto, la adición de tensoactivos en mezclas agua–alcano es de particular interés en el 
proceso de recuperación asistida [4]. 
Los tensoactivos, en inglés surfactants y su derivación en español surfactante, son agentes 
activos de superficie que disminuyen la tensión interfacial entre la solución de surfactantes y 
aceite residual. Estos agentes se adsorben en la superficie de la interface fluido–fluido. 
Comúnmente, la estructura de los tensoactivos se conforma de una cola no polar y una cabeza 
polar o iónica. Según las caracteristicas de sus partes hidrofóbicas e hidrofílicas, es el 
comportamiento del tensoactivo. Según la naturaleza de la cabeza, es la clasificación del 
tensoactivo, siendo: anionicos, catiónicos y no–iónicos o zwitterion [5]. 
Los tensoactivos de naturaleza no–iónica, a diferencia de los iónicos, no se ionizan en 
soluciones acuosas y tienen una doble función al actuar como agentes co–surfactantes y 
estabilizar las emulsiónes, además de ser más barata su elaboración en comparación con los 
surfactantes catiónicos [6]. 
Las mezclas agua–hidrocarburo, para fines de ingeniería y aplicación en diversas áreas, son 
de interés por ser sistemas compuestos por dos fases líquidas poco miscibles involucradas en 
una gran cantidadde procesos de separación o emulsificación. Este comportamiento se le 
atribuye a la naturaleza inherente a cada sustancia. Por un lado, las moléculas de hidrocarburo 
presentan interacciones por fuerzas de London, mientras que las moléculas de agua destacan 
por sus interacciones altamente direccionales y de corto alcance, los puentes de hidrógeno. 
Las parafinas o alcanos saturados como el n–hexano, n–heptano, n–octano, n–nonano y n–
decano son hidrocarburos sencillos de cadena larga que tienen una gran cantidad de 
aplicaciones en mezclas con otros hidrocarburos como en las gasolinas, e.g., nafta, diésel, 
queroseno, combustóleo, en la elaboración de aceites, lubricantes, como solventes e incluso 
11 
 
como materia prima o medio de reacción para determinadas sustancias. Los n–alcanos son 
sustancias orgánicas saturadas, compuestas por átomos de carbono e hidrógeno unidos por 
enlaces sencillos, poco solubles en agua y con un momento dipolar casi nulo, lo cual las hace 
fácilmente descriptibles como es uno de los primeros casos de estudio, el n–butano [7,8] 
(Figura 1.1). 
 
Figura 1.1 Representación de n–butano. 
Por otro lado, el agua es sin duda una de las más fascinantes; ocupando alrededor del 75 % 
de la superficie terrestre [9], siendo no solo responsable de sostener la vida, sino que también 
toma gran relevancia como solvente, soluto, reactante, catalizador, agente formador de 
estructuras en proteínas, ácidos nucleicos y células, etc. 
Muchos han sido los intentos por describir la tensión que actúa en la superficie de un líquido; 
por qué un líquido parece desafiar la gravedad dentro de un tubo, la peculiar curvatura que 
exhibe el menisco del mismo líquido en el capilar o la fuerza que actúa en algunos líquidos 
haciendo que al entrar en contacto se atraigan, son solo algunas preguntas que han motivado 
los estudios acerca del fenómeno de capilaridad. Los intentos para la descripción de este 
fenómeno remontan al siglo XVIII con modelos mecánicos [10,11]. Posteriormente, en el 
siglo XIX se dio paso a los métodos cuasi–termodinámicos que describen 
fenomenológicamente la naturaleza de los mismos [12,13]. 
Se define la tensión superficial como la fuerza que ejerce la superficie de un fluido con 
tendencia a minimizar el área de la superficie, debido a fuerzas intermoleculares, actuando 
como una fuerza en oposición al incremento del área del líquido. Van der Waals, en su seminal 
trabajo [1], define la tensión superficial como la energía capilar por unidad de área. Concluye 
que sus resultados coindicen con la definición propuesta por Gibbs, quien define la energía 
de exceso en la capa transitoria [13]. 
Van der Waals introdujo el término gradiente de la densidad de energía libre para describir la 
interfase líquido–vapor. Dicho término tomó relevancia en la transición de fases de la 
mecánica estadística desde los trabajos seminales por Ginzburg y Landau para 
superconductores [14]. Más adelante, Cahn y Hilliard [15] redescubrieron la teoría gradiente 
de van der Waals al aplicarla a mezclas binarias usando modelos de soluciones regulares. A 
la fecha, se han desarrollado diferentes formas de estimar o calcular la tensión superficial; 
entre ellas están el método de Paracoro [16], principio de estados correspondientes [4], teoría 
de funcional de las densidades clásico [17–20], teoría gradiente de interfases fluidas [4,21–
24], y simulaciones moleculares [25–27]. 
12 
 
En trabajos previos se ha determinado la tensión superficial e interfacial a través de la teoría 
gradiente de interfases fluidas acoplándola a diferentes ecuaciones de estado, e.g., ecuaciones 
tipo cúbicas y/o ecuaciones de la familia SAFT. Los sistemas de estudio al momento se han 
enfocado en n–alcanos, n–perfluoroalcanos, n–alcoholes y agua [28], así como mezclas de 
hidrocarburos con gases ligeros [29], agua–hidrocarburos y agua–aromáticos [23]. En 
algunos casos se ha logrado extender los cálculos para mezclas ternarias [30,31]. 
Por otro lado, las simulaciones moleculares han demostrado éxito para reproducir, predecir 
y “visualizar” fenómenos que se llevan a escala microscópica, difícilmente accesibles 
experimentalmente, siendo el tiempo de computo la única limitante para la obtención de 
resultados. Actualmente, con las computadoras que se cuenta hoy en día, es posible realizar 
estos cálculos de manera “rápida” inclusive en un ordenador personal. 
Una simulación molecular permite la cuantificación y visualización de características 
estructurales de las moléculas en la región no homogénea que es la interfase. La distribución 
de las moléculas a lo largo de la interfase se muestra a través de los perfiles de densidad. En 
mezclas, es posible visualizar la acumulación de algunas especies, como es el caso de los 
tensoactivos, tomando como referencia la posición en la interfase donde la concentración 
superficial de la especie de referencia es nula, denominada Superficie Divisoria de Gibbs 
(GDS). 
Experimentalmente, técnicas como la difracción y reflexión de rayos–X, reflexión de 
neutrones, microscopía de fuerza atómica y técnicas de espectroscopía óptica no–lineal como 
la espectroscopía vibracional de suma de frecuencias [32] son usadas para proveer 
información microscópica en interfases con una fase líquida presente. Sin embargo, estas 
técnicas no son usadas comúnmente en equilibrios líquido–vapor por detalles técnicos como 
el hecho de que las mediciones se realizan a condiciones de alto vacío o que se encuentren 
cubiertas. 
La información experimental de estudios realizados a la estructura de los alcanos en la 
interfase no es tan amplia (como en el caso del agua con su vapor [33,34] o sustancias con 
radicales OH en su estructura, como los alcoholes [32], incluso tensoactivos han sido sujeto 
de estudio en diversos trabajos [35]), viéndose reducidos los datos experimentales 
disponibles para alcanos de cadena larga (desde el hexadecano hasta el eicosano) [36]. 
Por otro lado, otros trabajos donde se han estudiado las orientaciones de las moléculas usando 
simulaciones moleculares han demostrado que existe un ordenamiento por parte de los 
alcanos en la interfase [37]. Este ordenamiento se intensifica a temperaturas bajas, y 
disminuye a temperaturas altas. 
Para este trabajo, el vector característico de los n–alcanos corresponde a un vector promedio 
a lo largo de la cadena, mientras que para el agua se usaron los vectores del momento dipolar 
del grupo OH y del vector normal al plano de la molécula HOH. 
La motivación para llevar a cabo este trabajo es la creciente necesidad de tecnologías para la 
emulsificación y desemulsificación en procesos industriales (extracción de petróleo asistida, 
13 
 
la remoción de contaminantes, extracción líquido–líquido) así como fenómenos biológicos 
(administración de fármacos). Se busca generar nuevos datos, debido a la escasa información 
dedicada a la inclusión de tensoactivos etoxilados a mezclas agua–hidrocarburos. 
 
1.1 Plan de tesis 
Este trabajo se divide de la siguiente forma: En el capítulo 2 se proporciona el marco teórico 
de las simulaciones moleculares usando la técnica de dinámica molecular, elementos de la 
teoría gradiente de fases fluidas y ecuaciones de estado para el cálculo de tensiones 
superficiales y las propiedades que caracterizan la interfase. En el capítulo 3 se detalla la 
metodología seguida para realizar las simulaciones, las funciones para describir a cada 
sustancia, la preparación de las cajas de simulación iniciales y las ecuaciones para el cálculo 
de propiedades. En el capítulo 4 se presentan los resultados y discusiones, producto de la 
interpretación de datos, para los hidrocarburos, agua, mezclas agua–hidrocarburo, agua–
tensoactivo y agua–hidrocarburo–tensoactivo. Finalmente, el capítulo 5 contiene las 
conclusiones y una serie de recomendaciones a considerar para futuros trabajos. 
 
14 
 
Referencias[1] Enhanced Oil Recovery, www.doe.gov. U.S. Department of Energy. 
[2] E. Tzimas, Enhanced Oil Recovery using Carbon Dioxide in the European Energy 
System, European Commission Joint Research Center. (2005) Retrieved 2012-11-01. 
[3] J.J. Sheng, Petroleum, 1 (2015) 97–105. 
[4] Y.X. Zuo, E.H. Stenby, Can. J. Chem. Eng., 75 (1977) 1130–1137. 
[5] S.B. Sandersen, PhD Thesis, Delft University of Technology, Delft, The Netherlands, 
2012. 
[6] C. Negin, S. Ali, Q. Xie, Petroleum, 3 (2017) 197–211. 
[7] J.P. Ryckaert, A. Bellemans, Chem. Phys. Lett., 30 (1975) 95–106. 
[8] J.P. Ryckaert, G. Ciccotti, J. Chem. Phys., 78 (1983) 7368–7374. 
[9] The Water Planet, www.nasa.gov. NASA. 
[10] H. Cavendish, Phil. Trans. Roy. Soc., 88 (1798). 
[11] F. Hauksbed, (1709) citado por Laplace (1806) (ref. 3, p. 3; Bowditch, ref. 3, p. 688). 
[12] J.D. van der Waals, J. Stat. Phys., 20 (1979) 200–244. 
[13] J.W. Gibbs, Trans. Connecticut Acad., 3 (1875–7) 108–248, 343–524. 
[14] L.D. Landau, Collected papers, (Oxford: Pergamon Press, 1965) 546–568. 
[15] J.W. Cahn, J.E. Hilliard, J. Chem. Phys., 28 (1958) 348–355. 
[16] D.B. MacLeaod, Trans. Faraday Soc., 9 (1923) 38–43. 
[17] R. Evans, Density functionals in the theory of nonuniform fluids, in: D. Henderson (Ed.), 
Fundamentals of Inhomogeneous Fluids, Dekker, New York, 1992. 
[18] G.J. Gloor, G. Jackson, F.J. Blas, E. Martín del Río, E. de Miguel, J. Chem. Phys., 121 
(2004) 12740–12759. 
[19] Z. Li, A. Firoozabadi, J. Chem. Phys., 130 (2009) 154108. 
[20] F. Llovell, A. Galindo, F.J. Blas, G. Jackson, J. Chem. Phys., 133 (2010) 24704. 
[21] V. Bongiorno, L.E. Scriven, H.T. Davis, J. Colloid Interface Sci., 57 (1976) 462–475. 
[22] C.I. Poser, I.C. Sanchez, Macromolecules, 14 (1981) 361–370. 
[23] M. Medeiros, Fluid Phase Equilib., 392 (2015) 65–73. 
[24] J.M. Garrido, A. Mejía, M.M. Piñeiro, F.J. Blass, E.A. Müller, Am. Inst. Chem. Eng. J., 
62 (2016) 1781–1794. 
[25] E.A. Müller, A. Mejía, Fluid Phase Equilib., 282 (2009) 68–81. 
[26] C. Miqueu, J.M. Míguez, M.M. Piñeiro, T. Lafitte, B. Mendiboure, J. Phys. Chem. B, 
115 (2011) 9618–9625. 
[27] A. Mejía, C. Herdes, E.A. Müller, Ind. Eng. Chem. Res., 53 (2014) 4131–4141. 
[28] M.B Oliveira, I.M. Marrucho, J.A.P. Coutinho, A.J. Quemeida, Fluid Phase Equilib., 
267 (2008) 83–91. 
[29] S.T Liu, D. Fu, J.Y. Lu, Ind. Eng. Chem. Res., 48 (2009) 10734–10739. 
[30] C. Miqueu, B. Mendiboure, C. Graciaa, J. Lachaise, Fluid Phase Equilib., 218 (2004) 
189–203. 
[31] E. Schäfer, F. Horbach, S. Enders, J. Chem. Eng. Data, 59 (2014) 3003–3016. 
[32] P.B. Miranda, Y.R. Shen, J. Phys. Chem. B, 103 (1999) 3292–3307. 
[33] Q. Du, R. Superfine, E. Freysz, and Y. R. Shen. Phys. Rev. Lett., 70 (1993) 2313–2316. 
http://www.doe.gov/
http://www.nasa.gov/
15 
 
[34] Y. Fan, X. Chen, L. Yang, P.S. Cremer, Y.Q. Gao. J. Phys. Chem. B, 113 (2009) 11672–
11679. 
[35] E. Tyrode, C.M. Johnson, A. Kumpulainen, M.W. Rutland, P.M. Claesson, J. Am. Chem. 
Soc., 127 (2005) 16848–16859. 
[36] G.A. Sefler, Q. Du, P.B. Miranda, Y.R. Shen, Chem. Phys. Lett., 235 (1995) 347–354. 
[37] P. Smith, R.M. Lynden–Bell, J.C. Earnshaw, W. Smith, Mol. Phys., 96 (1999) 249–257. 
 
16 
 
Objetivos 
El objetivo principal es calcular las tensiones interfaciales de mezclas agua–hidrocarburo (n–
hexano, n–heptano, n–octano, n–nonano y n–decano) en presencia de tensoactivos 
pertenecientes a la familia de los alcoholes etoxilados (C12E2, C12E3, C12E4 y C12E5) y 
cuantificar el efecto de la concentración y el tamaño de la cabeza hidrofílica de los diferentes 
tensoactivos sobre la tensión interfacial, distribución, orientación y movilidad de las 
moléculas de agua e hidrocarburos en la región interfacial. 
Se validaron los datos generados provenientes de los modelos TraPPE–UA, OPLS–UA, 
TIP4P y TIP4P/2005; comparando con mediciones experimentales de tensión superficial para 
cada sustancia pura, i.e., n–hexano, n–heptano, n–octano, n–nonano, n–decano y agua, así 
como la tensión interfacial de las mezclas binarias agua–hexano, agua–heptano, agua–
octano, agua–nonano, agua–decano, agua–C12E2, agua–C12E3, agua–C12E4 y agua–C12E5. 
Además de las tensiones superficiales e interfaciales, se desea caracterizar la estructura 
microscópica a través de los perfiles de densidad local, conformación y orientación molecular 
en la región interfacial a través de la distancia de extremo a extremo, radio de giro y 
parámetro de orden orientacional en la región interfacial de las sustancias puras y mezclas. 
Al final, se pretender establecer una relación entre las propiedades superficiales, i.e., tensión 
superficial e interfacial, y las propiedades microscópicas. 
Se caracterizaron las propiedades dinámicas de las especies en las mezclas, i.e., coeficiente 
de difusión lateral de las sustancias puras y mezclas, así como cambios energéticos en la 
interfase usando las simulaciones moleculares como herramienta de interpretación. 
 
 
17 
 
Hipótesis 
A partir de simulaciones moleculares y cálculos provenientes de la teoría gradiente de 
interfases fluidas, es posible interpretar los fenómenos superficiales que se llevan a cabo en 
la región interfacial de fases fluidas de especies moleculares simples (n–alcanos) y complejas 
(agua) a nivel microscópico y su efecto sobre la tensión superficial e interfacial. 
La interfase es una región especial con un espesor, en donde se llevan a cabo, principalmente, 
cambios orientacionales. Es posible comparar los espesores, distribución y tensión superficial 
predicha por la teoría gradiente de interfases fluidas con la dinámica molecular. 
Los cambios energéticos entre fases fluidas en contacto, conformaciones y orientaciones 
moleculares, así como las propiedades dinámicas de cada especie son funciones locales que 
dependen de la distribución y orientación molecular de cada especie, según su posición en la 
interfase del sistema. 
 
18 
 
2 Marco teórico 
El objetivo de usar simulaciones de dinámica molecular (DM) es obtener las propiedades 
mecánicas promediadas de un sistema dinámico, compuesto por N partículas que siguen las 
leyes de la física, de manera numérica en un periodo de tiempo suficientemente largo. A partir 
de la información generada a nivel microscópico se obtienen propiedades a nivel 
macroscópico (e.g., temperatura, presión, energía interna, etc.), tarea de la mecánica 
estadística. 
Como sugiere Gibbs [1], el estado dinámico de un sistema se puede especificar detallando la 
información mecánica del sistema en un determinado punto en un espacio de 2n dimensiones, 
i.e., posición y momento de un sistema con n grados de libertad. A esta construcción 
matemática se le denomina espacio fásico. 
El valor instantáneo de cualquier propiedad termodinámica (e.g., temperatura, presión, 
volumen, energía interna), G, está determinado por el estado dinámico del sistema, i.e., 
G(p,q), donde p y q son las coordenadas de momento y posición generalizadas para cada 
partícula que conforma el sistema de estudio. Debido a que el sistema evoluciona con el 
tiempo, p y q lo hacen al igual que G(p,q), por lo tanto es razonable asumir que la propiedad 
macroscópica observable G(p,q) en realidad es un promedio temporal en un intervalo amplio 
de tiempo [2]. 
obs
obs
obs temp temp
0
obs
( ( ), ( ))
1
lim ( ( ), ( )) .
t
t
G G G p t q t
G p t q t dt
t
 
 
 (2.1) 
Para obtener la precisión deseada no se puede esperar que la Ec. 2.1 se resuelva para tiempos 
infinitos, por lo que los promedios son por intervalos de tiempo finito. La DM resuelve las 
ecuaciones de movimiento por intervalos de tiempo, paso a paso, de manera que 
obs
obs temp
1obs
1
( ( ), ( )) ;G G G p t q t

 
   (2.2) 
donde  es el índice de la suma correspondiente al número de intervalos temporales. 
Actualmente las computadoras son capaces de llevar a cabo tal cantidad de cálculos, y es 
factible explorar de manera satisfactoria el espacio fásico delsistema para obtener promedios 
temporales en un tiempo de computo razonable, así como lograr igualdad de condiciones 
entre simulaciones con parámetros macroscópicos y diferentes condiciones iniciales. 
Debido a la complejidad que representa el cálculo de las propiedades G(p,q) para un gran 
número de moléculas, Gibbs sugirió reemplazar el promedio temporal por el promedio de 
colectivo [3]. El colectivo es una colección de puntos p y q en el espacio fásico distribuidos 
de acuerdo a una densidad de probabilidad (p,q). Esta función se determina al fijar un 
juego de parámetros macroscópicos (NPT, NVT, etc.). 
19 
 
Las ecuaciones de movimiento de Newton generan sucesiones de estados conforme a la 
función de distribución del colectivo microcanónico, NVE. Para otros colectivos se requiere 
que las ecuaciones de movimiento cumplan con las siguientes condiciones: 
i) La densidad de probabilidad del colectivo no debe cambiar con el tiempo. 
ii) Cualquier distribución inicial tiende a una solución estacionaria en el transcurso de la 
simulación. 
iii) Se debe sostener la hipótesis ergodica, aunque no se pueda probar para sistemas 
realistas. 
De esta manera, cualquier estado inicial razonable es capaz de generar una sucesión de 
puntos, a largo plazo, de acuerdo a la densidad de probabilidad deseada. 
 
2.1 Dinámica molecular 
Las bases teóricas detrás de la dinámica molecular (DM) involucran trabajos producidos por 
grandes personajes que han contribuido ampliamente en el campo de la mecánica analítica, 
tales como Euler, Hamilton, Lagrange, Newton [1]. Algunos de los resultados contienen 
observaciones fundamentales acerca de la naturaleza física, mientras que otros son 
reformulaciones elegantes que invocan nuevos desarrollos teóricos. 
La DM es un método determinista que requiere el conocimiento de la posición y momento 
de una partícula para determinar su pasado, presente y futuro (trayectoria), lo cual requiere 
la segunda ley de Newton. En caso de que las moléculas sean rígidas es necesario el uso de 
las ecuaciones de Euler, posiblemente expresadas en forma de los cuaterniones, propuestos 
por Hamilton. Moléculas con grados de libertad internos sujetas a restricciones estructurales 
requieren incorporar restricciones geométricas, usando el método de Lagrange, a las 
ecuaciones de movimiento. Por otro lado, las ecuaciones de movimiento de Newton (EdM), 
en la DM, corresponden al colectivo microcanónico de la termodinámica estadística; sin 
embargo, si se requiere que las condiciones de simulación igualen las experimentales (i.e., 
temperatura y/o presión constante) se deben modificar las EdM para reproducir las variables 
macroscópicas, aunque la trayectoria ya no representa la solución a las ecuaciones de 
Newton. 
Debido a la complejidad que representa la solución de las EdM y el número de partículas que 
conforman el sistema, las EdM se resuelven numéricamente, por lo que se debe evaluar 
cuidadosamente la metodología para resolver las ecuaciones, ya que incorpora un error 
inherente al método numérico. 
Debido a que algunas interacciones intermoleculares son sensibles a perturbaciones, cambios 
pequeños pueden traer efectos que se magnifican de manera exponencial, lo cual limita los 
periodos de simulación. Todos estos puntos son clave para calcular trayectorias con precisión 
moderada, que preserven la energía del sistema, que las EdM sean reversibles en tiempo [5]. 
20 
 
2.1.1 Ecuaciones de movimiento 
La solución de las ecuaciones de movimiento es el corazón de las simulaciones de DM 
clásica, i.e., la integración numérica de las ecuaciones de movimiento de Newton para un 
sistema compuesto de N–partículas interactuantes: 
2
2
 ,ii im
t



r
F (2.3) 
donde mi es la masa de la partícula i, ri el vector de posición y Fi la fuerza que actúa en la 
partícula i. La fuerza está dada por 
( )
i
N
i U rF r (2.4) 
donde U es la función potencial que representa la energía potencial del sistema. La función 
potencial, en su forma más sencilla, se puede aproximar como una suma de interacciones por 
pares: 
1
( )
N N
i j
ij
i
U u r
 
 (2.5) 
donde u es el potencial por pares i–j, rij es la distancia de separación entre dos partículas 
definida como rij=|ri – rj|. Existen diversos métodos de integración numérica que resuelven 
las ecuaciones de movimiento paso a paso; sin embargo, debido a la cantidad de partículas 
dentro de la caja de simulación, el tiempo de simulación requerido para el muestreo de las 
propiedades y la complejidad del mismo método, se opta por usar métodos sencillos que sean 
capaces de conservar la energía y capaces de reproducir ciertas correlaciones de tiempo y 
espacio con un determinado grado de precisión [6], siendo el método más común el algoritmo 
salto de rana que se describe a continuación. 
2.1.2 Algoritmo salto de rana 
Los algoritmos de salto de rana [5] son métodos de integración numérica que evalúan las 
velocidades a la mitad del intervalo de tiempo propuesto y con estas evalúan la nueva 
posición. Este método tiene la ventaja de ser un método de fácil programación, los 
requerimientos de almacenamiento son mínimos, tienen una precisión del orden de magnitud 
~O(t3) y preservan mejor cualquier cantidad conservada (energía, momento angular) al ser 
reversibles. A continuación, se describe el algoritmo del método de Verlet que se puede 
reescribir de forma intuitiva para obtener el método de salgo de rana. La serie de Taylor de 
la variable de coordenada x(t) truncada en el tercer término es 
2
3( ) ( ) ( ) ( ) ( ) ;
2
h
x t h x t hx t x t O h     (2.6) 
donde t es el tiempo actual y h≡t es el intervalo de tiempo transcurrido. La notación punto 
se refiere al componente velocidad y aceleración (ẋ y ẍ, respectivamente) en función de t, 
aunque, ẍ es una función conocida que depende de las coordenadas en el tiempo t. Sumando 
la expresión que corresponde a un paso inverso x(t – h) a la Ec. 2.6, obtenemos: 
21 
 
2 4( ) 2 ( ) ( ) ( ) ( )x t h x t x t h h x t O h      (2.7) 
Lo cual cambia el error de truncamiento a O(h4) por la cancelación del término h3. La 
velocidad se obtiene usando el método de punto–medio. 
2( ) ( )( ) ( )
2
x t h x t h
x t O h
h
  
  (2.8) 
Reescribiendo la Ec. 2.6: 
3( ) ( ) ( ) ( ) ( )
2
h
x t h x t h x t x t O h
 
     
 
 (2.9) 
( / 2) ( ) ( )
2
h
x t h x t x t   (2.10) 
Restando a la Ec. 2.10 la velocidad evaluada en t – h/2 se obtiene 
( / 2) ( / 2) ( )x t h x t h hx t    (2.11) 
Sustituyendo la Ec. 2.10 en la Ec. 2.9 
( ) ( ) ( / 2)x t h x t hx t h    (2.12) 
Que representan las ecuaciones del algoritmo de salto de rana. Reescribiendo las ecuaciones 
correspondientes para la posición y velocidad en el instante t en forma vectorial: 
1 1
( )
2 2
i i it t t t t t
   
         
   
v v a (2.13) 
1
 ( ;
2
) ( )i i it t tt t t
 
    

 

 vr r (2.14) 
donde r es el vector posición, v es el vector velocidad y a el vector aceleración de la partícula 
i, con un intervalo de tiempo equivalente a la mitad del intervalo de tiempo t. Este método 
requiere calcular la velocidad en el tiempo actual 
1 1 1
( )
2 2 2
    
         
    
i i it t t t tv v v (2.15) 
Esta velocidad es necesaria para el cálculo de la energía total del sistema. Una vez evaluadas 
las propiedades de interés la posición actual se usa para calcular la nueva posición de la 
partícula al tiempo t+t. 
2.1.3 Condiciones de frontera 
Establecer el tamaño de la caja de simulación es una tarea difícil, ya que un sistema grande 
producirá resultados más fieles a sistemas de tamaño infinito; sin embargo, esto traerá como 
consecuencia que el tiempo de simulación sea mucho mayor. Por otro lado, el fijar paredes 
rígidas enla caja de simulación para evitar que las partículas salgan, provoca una 
perturbación en las proximidades a la frontera, cambiando el comportamiento en la cercanía 
22 
 
del bulto de la fase fluida. A menos que se desee estudiar la interacción de partículas con una 
interfase sólida es mejor eliminar las paredes. 
Las condiciones de frontera periódicas (CFP) es una forma de eliminar las barreras físicas 
impuestas a un sistema y considerar que la caja de simulación tiene un tamaño “infinito”, al 
realizar copias idénticas en todas direcciones de la caja de simulación (Figura 2.1). Las 
consecuencias por implementar las CFP son: 
i) Al salir una partícula de la caja de simulación, una copia entra por el lado opuesto de 
la caja. 
ii) Una partícula es capaz de interactuar con partículas dentro de la caja de simulación y 
sus respectivas copias dentro de una esfera virtual con un radio, radio de corte rc. 
 
 
Figura 2.1 Representación esquemática de condiciones de frontera periódica con partículas 
dentro de la caja de simulación (círculo rojo), interactuantes dentro del radio de corte rc 
(círculo lleno) con copias de partículas (círculo con trama). 
Durante la evolución temporal del sistema, es necesario considerar el efecto de las CFP 
después de cada paso de integración, por lo que se debe verificar si las partículas salieron de 
la caja de simulación y de ser así ingresarlas por el extremo en dirección opuesta. De igual 
manera, al realizar el análisis de la trayectoria de la simulación, se debe considerar el efecto 
de las CFP. 
Para considerar las interacciones de las partículas con otras que se encuentran fuera del radio 
de corte es necesario realizar correcciones a las interacciones de largo alcance a la energía y 
presión. La contribución por dichas correcciones a la energía y presión del sistema son 
2
la 2 ( )
cr
U u r drN r 

  (2.16) 
2
la 2 ( )
cr
N r r dr  

   (2.17) 
23 
 
donde Ula es la corrección a la energía potencial de largo alcance,  es el número de densidad 
promedio y la es la corrección al término del virial (ecuación que relaciona la energía 
cinética con la energía potencial) de largo alcance [7]. 
2.1.4 Interacciones Coulombicas (Sumas de Ewald) 
Un sistema compuesto por N cargas localizadas en r1, r2, r3, …, rN su energía potencial 
Coulómbica es 
0
1
4
i j
iji j
q q
U
r
 (2.18) 
donde qi es la carga de la partícula i. Suponiendo que las cargas están sujetas a condiciones 
de frontera periódica, CFP, descritas por vectores repetidos (c1, c2, c3) conformando una celda 
unitaria de mayores dimensions (supercelda), las cargas iónicas qi también están repetidas 
ubicadas en ri+n1c1+n2c2+n3c3, donde n1, n2 y n3 son números enteros arbitrarios. Entonces 
la energía potencial de la supercelda es 
0
1
'
4
i j
iji j
q q
U
r L


 
n n
 (2.19) 
donde para una supercelda cúbica L es la longitud de la supercelda con L=|c1|=|c2|=|c3| y los 
vectores n representan la imagen periódica. La suma excluye al termino j=i si y solo si n=0. 
La Ec. 2.19 es una suma que converge lentamente y es condicionalmente convergente [8], lo 
cual implica que depende del orden de la suma. 
2.1.4.1 Perspectiva física 
El método de Ewald [8,9] evalúa la suma al transformarla en una suma que converge 
rápidamente y de manera absoluta, lo cual implica el separar la distribución de cargas en tres 
términos. En este apartado se resume en unos comentarios los puntos importantes de la 
técnica. 
La clave para calcular de manera rápida las energías y fuerzas consiste en separar la suma de 
la siguiente forma 
1 ( ) 1 ( )f r f r
r r r

  (2.20) 
donde f(r) es una función que decae rápidamente y corresponde a la función error 
complementaria. La idea básica consiste en separar la parte que varía rápidamente para 
valores pequeños de r de manera que decaigan rápido y su contribución más allá del radio de 
corte rc sea despreciable y la parte suave para valores grandes de r, la cual se manipula en el 
espacio de Fourier [10]. Ewald sugiere elegir un factor de convergencia g(s,n)=e-s|n‧n|, tal que 
la suma se convierte en 
2
1
2
2 2
0
4
,
2
erfc( )1
'
2
1 4
( )
2 k
k i j ij
i
i j iij
q q r L
q
r
U k e
V L

 



  
  
 
 

 
nk
n
n
 (2.21) 
24 
 
donde  es denominado el parámetro de separación de Ewald que separa la parte real y 
recíproca. La elección de qué tan grande se desea la supercelda y la rapidez con que decrece 
la contribución de los términos de la suma es de acuerdo con la exactitud que se desee 
alcanzar. 
Existen variaciones de este método para mejorar el desempeño de la suma reciproca, siendo 
uno de los más populares el método de malla de partículas de Ewald [11], en inglés PME, el 
cual asigna las cargas a una cuadrícula, interpolando con un spline. La cuadrícula es 
transformada al espacio Fourier y la energía reciproca se obtiene a través de una suma en la 
cuadricula del k–espacio. El potencial en cada punto de la cuadrícula se obtiene 
transformando al espacio real y usando los factores de interpolación se obtiene la fuerza en 
cada átomo [11]. 
2.1.5 Funciones potenciales 
Comúnmente, la energía potencial de un grupo de partículas se divide en dos contribuciones: 
externas que se originan por interacciones con partículas pertenecientes a otras moléculas y 
campos externos e internas que son el resultado de interacciones entre partículas que 
pertenecen a la misma molécula. 
La energía externa se divide en cuatro contribuciones: 
ext el pol disp repU U U U U    (2.22) 
donde Uel es la energía potencial electrostática, Upol la energía de polarización, Udisp es la 
energía de dispersión y Urep la energía de repulsión. Por otro lado, la energía interna 
(intramolecular) se expresa como una suma de contribuciones 
intra enlace angular diedro no enlace ,U U U U U     (2.23) 
donde Uenlace es la energía asociada a las variaciones de longitud de enlace, Uangular es la 
energía asociada a la variación del ángulo formado por dos enlaces sucesivos, Udiedro es la 
energía causada por la variación del ángulo diedro formado por cuatro partículas sucesivas 
en una cadena y Uno–enlace es la energía resultante entre partículas pertenecientes a la misma 
molécula y separadas por más de tres enlaces. En la Figura 2.2 se muestra la interacción de 
no–enlace entre los sitios 1–6 de una cadena. 
 
Figura 2.2 Representación de cadena con interacciones 1–6. 
El campo de fuerza (F), resultado del gradiente de la función potencial, es una representación 
del campo vectorial que actúa en una partícula por interactuar con otros cuerpos o campos 
externos. En la actualidad, existe un gran número de potenciales disponibles para el cálculo 
de propiedad [12–14]; sin embargo, la selección del potencial más adecuado no es tarea 
1 
6 
25 
 
sencilla, ya que ninguna función es capaz de predecir todas las propiedades termodinámicas 
de forma cuantitativa, motivo por el cual existen tantos campos de fuerza, por lo que se debe 
acotar las propiedades de interés. Dentro de los puntos deseables en una función potencial 
están: 
• Que haya reproducido previamente las propiedades de interés exitosamente. 
• Que sea una función potencial con parámetros transferibles a sustancias diferentes a 
las reportadas. 
• Que se pueda ajustar fácilmente para reproducir otras propiedades. 
• Que sea compatible con otras funciones potenciales en mezclas. 
• Que las mezclas reproduzcan y/o predigan propiedades termodinámicas. 
Al representar moléculas poliatómicas, dos enfoques se pueden usar para representar las 
fuerzas de dispersión y repulsión: 
• Se asigna un centro de interacción a cada átomo de la molécula ubicado en su núcleo. 
Dicha metodología se le denomina “Todos los Átomos”, en ingles AA [15]. 
• Se asigna un centro de interacción a un grupode átomos como los grupos CH, CH2 o 
CH3. A esta metodología se le denomina “Átomos Unidos”, en ingles UA [16]. 
La diferencia entre una y otra radica en la aplicación que se da a cada metodología. Por un 
lado, el método AA brinda información explícita de la geometría molecular y su estructura, 
pero requiere un mayor tiempo de computo. El método UA reduce el tiempo de cómputo al 
considerar a un grupo de átomos como un pseudo–átomo. Aunque se pierde detalle 
molecular, el método UA ha reportado [17] resultados exitosos para equilibrios de fases 
fluidas y problemas de adsorción. 
2.1.5.1 Algoritmos de restricción 
Los algoritmos de restricción son métodos para satisfacer las ecuaciones de movimiento de 
Newton de cuerpos rígidos compuestos por masas puntuales [5]. Estos algoritmos se 
aseguran de mantener una distancia fija entre las masas. Los pasos que involucran consisten 
en: 
i) Elegir coordenadas sin restricciones. 
ii) Introducir fuerzas explicitas restrictivas. 
iii) Minimizar las fuerzas restrictivas utilizando la técnica de multiplicadores de 
Lagrange o proyecciones. 
Diferentes algoritmos para cumplir con esta tarea existen; entre ellos se encuentran los 
algoritmos SHAKE [18], RATTLE [19], SETTLE [20] y LINCS [21]. 
2.1.5.2 Potenciales del agua 
Existe un gran número de funciones potenciales usadas para determinar las propiedades del 
agua, e.g., SPC/E [22], TIP4P [23], TIP4P/2005 [24] y TIP4P/Ew [25], por mencionar 
algunos, los cuales suponen una estructura rígida del agua. El modelo SPC/E considera 3 
sitios de interacción centrados en los átomos de hidrógeno y oxígeno, mientras que los 
modelos pertenecientes a la familia TIP4P son modelos que consideran 4 sitios de 
26 
 
interacción: dos átomos de hidrógeno, un átomo de oxígeno y un sitio virtual que únicamente 
representa una carga negativa (Figura 2.3). 
 
 
Figura 2.3 Representación de la molécula de agua con 3 (panel a) y 4 sitios de interacción 
(panel b). Los átomos están representados por esferas de colores siendo: átomo de hidrógeno 
(gris), átomo de oxígeno (rojo), sitio virtual (amarillo). 
El potencial de los modelos se basa en los potenciales de Coulomb y Lennard–Jones: 
12 6
no enlac
0
e
1
4 ;
4
i j
iji j ij ij
q q
U
r r r
 




    
              
 (2.24) 
donde 0 es la permitividad del vacío, qi es la carga de la especie i, rij es la distancia entre las 
partículas i y j,  es la distancia en la que el potencial entre partículas es 0 y  es la 
profundidad del potencial. Los parámetros de los potenciales se muestran en la Tabla 2.1. 
Tabla 2.1 Parámetros de funciones potencial para molécula de agua. 
 SPC/E [22] TIP4P [23] TIP4P/2005 [24] TIP4P/Ew [25] 
, nm 0.3166 0.3154 0.3159 0.3164 
/kB, K 78.18 78.00 93.20 81.90 
, nm – – – – 
/kB, K – – – – 
dOH, nm 0.10000 0.09572 0.09572 0.09572 
dOM, nm – 0.01500 0.01546 0.0125 
qO/e
– +0.4238 +0.5200 +0.5564 +0.52422 
 –0.8476 – – – 
qM/e
– – –1.04000 –1.11280 –1.04844 
, ° 109.47 104.52 104.52 104.52 
, ° – 52.26 52.26 52.26 
 
donde kB es la constante de Boltzmann, dOH es la distancia entre el átomo de oxígeno y el 
átomo de hidrógeno, dOM es la distancia entre el átomo de oxígeno y el sitio virtual M,  es 
el ángulo entre las especies HOH y  el ángulo entre las especies MOH. Todas las moléculas 
de agua mantuvieron su rigidez con el algoritmo SETTLE [20]. 
a b 
27 
 
2.1.5.3 Potenciales OPLS–UA y TraPPE–UA 
La función potencial “Potenciales Optimizados para Simulación de Líquidos”, en inglés 
OPLS, propuesta por Jorgensen et al. [16], pertenece a una familia de funciones para la 
simulación de sustancias de naturaleza orgánica y sistemas bioquímicos en fase líquida, 
parametrizadas para reproducir propiedades termodinámicas y estructurales en el bulto de 
fases. Dentro de la familia de potenciales OPLS se pueden representar los sitios de 
interacción usando los enfoques OPLS–UA y OPLS–AA, siendo el potencial OPLS–UA el 
usado en este trabajo. 
Por otro lado, la función potencial “Potenciales Transferibles para Equilibrios de Fase”, en 
inglés TraPPE–UA, propuesta por Martin y Siepmann [26] se usa para predecir las 
propiedades de hidrocarburos, basada en la función potencial OPLS–UA, logrando reducir 
el tiempo de computo al considerar interacciones tipo átomo unido. A diferencia de funciones 
potenciales similares (OPLS–UA y SKS) el potencial TraPPE–UA se ajustó a curvas de 
coexistencia de equilibrio. Müller y Mejía [17] mostraron que el campo de fuerza TraPPE–
UA logra reproducir densidades experimentales y tensiones superficiales de hidrocarburos 
lineales con desviación baja. 
La forma funcional de los potenciales OPLS–UA y TraPPE–UA es 
  no enlace ángulo– diedro ,NU U U U  r (2.25) 
donde no se considera un término referente a la energía de enlace por los enlaces rígidos de 
la molécula. A continuación, se presentan las interacciones de no–enlace para el potencial 
aditivo por pares de Coulomb y Lennard–Jones 12–6 (Ec. 2.24). Los parámetros de los 
pseudo–átomos se muestran en la Tabla 2.2 
Tabla 2.2 Parámetros de no–enlace para los potenciales OPLS–UA [16,27] y TraPPE–UA 
[26]. 
Potencial Pseudo–átomo /kB, K , nm q/e
– 
OPLS–UA 
CH2 (éter) 59.38 0.3800 +0.250 
CH2 (hidroxilo) 59.38 0.3905 +0.265 
O (éter) 85.55 0.3000 –0.500 
O (hidroxilo) 85.55 0.3070 –0.700 
H (hidroxilo) 0.00 0.0000 +0.435 
TraPPE–UA 
CH3 (n–alcano) 98.00 0.3750 0.0 
CH2 (n–alcano) 46.00 0.3950 0.0 
 
donde e– es la carga fundamental del electrón. Los parámetros de interacción de no–enlace 
entre especies diferentes se estimaron usando las reglas de combinación de Lorentz–
Berthelot que consideran un promedio geométrico para ij y un promedio aritmético para ij. 
A continuación, se presenta el potencial armónico angular formado por tres pseudo–átomos 
(Figura 2.4) 
28 
 
 ángu 0lo
2
2
k
U    (2.26) 
 
Figura 2.4 Representación esquemática del potencial armónico angular. 
Los parámetros entre diferentes especies para el potencial armónico angular se muestran en 
la Tabla 2.3. 
Tabla 2.3 Parámetros del potencial de armónico angular para los potenciales OPLS–UA 
[16,27] y TraPPE–UA [26]. 
∠ 0, ° k/kB, K rad–2 
CHx–CHx–CHx 114.0 62500 
CHx–CHx–O
éter 112.0 50300 
CHx–O
éter–CHx 112.0 60400 
Oéter–CHx–CHx 112.0 62500 
CHx–CHx–O
OH 108.0 50402.8 
CHx–O
OH–HOH 108.5 55403.1 
 
donde ∠ es el ángulo formado entre las especies listadas, k es la constante del potencial,  
es el ángulo formado por las especies y 0 es el ángulo de equilibrio. 
El potencial de diedro representa el ángulo entre dos planos que intersectan (Figura 2.5). 
 
Figura 2.5 Representación esquemática del potencial de diedro. 
El potencial de diedro se puede expresar en forma de una serie de Fourier (Ec. 2.27). Las 
constantes usadas en este trabajo se muestran en la Tabla 2.4. 
     rot.int. 1 32
1
1 cos( ) 1 cos(2 ) 1 cos(3 )
2
( )U c c c        (2.27) 
29 
 
   
12 6
1
diedro rot.int.
0
,4
,
1
4
4
i
i j
i j
i i j ijj iji
q q
U U f
r r
r
r
 
 



                        
  (2.28) 
Tabla 2.4 Parámetros del potencial de diedro para los potenciales OPLS–UA [16,27] y 
TraPPE–UA [26]. 
∠ c1/kB, K c2/kB, K c3/kB, K 
CHx–CHx–CHx–CHx 355.03 –68.19 791.32 
CHx–CHx–CHx–O 176.64 –53.34 769.97 
CHx–CHx–O–CHx 725.36 –163.73 558.18 
CHx–CHx–O–H 209.85 –29.19 187.96 
 
donde  es el ángulo del diedro, c1, c2 y c3 son constantes del potencial. Las partículas i–j 
separadas por más de tres enlaces presentan interacciones de no–enlace, similar a la Ec. 2.24
. A estas interacciones se les denomina interacciones 1–4 y normalmente son escaladas por 
un factor fij. En el caso de las partículas descritas por el modelo TraPPE–UA las constantes 
de la Tabla 2.4 están parametrizadas para incluir estasinteracciones, i.e., fij=0. Por otro lado, 
Aunque el modelo OPLS–UA está parametrizado de forma similar, este considera 
únicamente las interacciones 1–5 para alcoholes y éteres (sustancias base para la cabeza 
hidrofílica). La Figura 2.6 representa las interacciones 1–5 de no–enlace en la molécula 
etilénglicol. 
 
 
Figura 2.6 Sitios interactuantes 1–5 entre las especies H–O (a) y O–CH3 (b) en el 
dietilenglicol. 
Las interacciones 1–5 correspondientes los parámetros de no–enlace se muestran en la Tabla 
2.5. 
Tabla 2.5 Parámetros de no–enlace por interacción entre pseudo–átomos 1–5 [16,27]. 
Par 1–5 Segmento /kB, K , nm 
C–C R–CH2–O
éter–CH2–R 4.026 0.2745 
C–O R–CH2–O
éter–CH2–COH 4.026 0.2745 
O–H ROéter–CH2–CH2–O–H 4.026 0.2600 
a 
a’ 
b’ 
b 
a 
b b’ a’ 
30 
 
 
Las moléculas representadas por los potenciales OPLS–UA y TraPPE–UA mantienen una 
distancia fija entre partículas con enlaces rígidos usando el algoritmo LINCS [21]. En la 
Tabla 2.6 se muestran las distancias entre cada especie. 
Tabla 2.6 Distancias intermoleculares entre átomos y/o pseudo–átomos OPLS–UA [16,27] y 
TraPPE–UA [26]. 
Pares de pseudo–átomos d / nm 
CHx–CHx 0.1540 
CH2–O
éter 0.1410 
CH2–O
OH 0.1430 
O–H 0.0945 
 
 
2.2 Propiedades interfaciales 
 
La tensión superficial e interfacial (energía libre de exceso) es la propiedad que se puede 
comparar directamente con datos experimentales. La energía libre de exceso, al igual que las 
demás propiedades de exceso, es una diferencia entre las propiedades de todo el sistema y 
las propiedades del bulto de cada fase sin que estas manifiesten algún cambio hasta la 
superficie divisoria [28]. 
El parámetro de orden mide el grado de ordenamiento a través de las fronteras en un sistema 
con transición de fases. Según la transición de fases es la definición del parámetro de orden. 
En este trabajo se usó un criterio similar al de cristales líquidos donde la orientación 
molecular respecto a un eje de laboratorio se cuantifica (parámetro de orden orientacional). 
La distancia de extremo a extremo de una cadena y el radio de giro proveen información 
acerca del tamaño de cadenas; sin embargo, en este trabajo se usaron para verificar si existe 
algún cambio conformacional en la molécula por efectos de la transición de fase y verificar 
la orientación preferencial de la molécula. 
La dinámica de las moléculas se puede seguir vía el coeficiente de difusión. Para sistemas 
no–homogéneos el coeficiente de difusión cambia según la posición de las partículas, en este 
caso su posición en la interfase. Para cuantificar la movilidad a lo largo de la interfase se 
mide la difusión en el plano normal a la interfase para cada segmento (difusión paralela o 
lateral). 
El cambio de energía libre cuantifica la diferencia energética de un estado a otro. En este 
trabajo, se determinó el cambio energético de una molécula de n–alcano por migrar de una 
fase líquida a vapor. Usando el muestreo de sombrilla, es posible seguir cómo cambia la 
diferencia de energía libre según su posición en la interfase [29]. Existe una relación entre el 
31 
 
cambio de energía libre y la energía de exceso al referirse ambas a cambios energéticos en la 
región interfacial. 
2.2.1 Perfiles de densidad 
Para calcular los perfiles de densidad locales, se dividió la caja de simulación en segmentos 
con espesor de 0.1 nm en la dirección del eje z. La densidad promedio en cada sección 
transversal de la caja se determinó tomando como referencia la posición del centro de masa 
de cada molécula. 
Existen diferentes definiciones para describir los perfiles de densidad en las simulaciones 
moleculares; sin embargo, para los equilibrios líquido–vapor es común encontrar en la 
literatura que estos perfiles son mejor descritos por funciones hiperbólicas [30]. 
    0L V L V
1 1
( ) tanh
2 2
z z
z
d
    
 
     
 
 (2.29) 
donde V y L son la densidad de vapor y líquido en el bulto de la fase, respectivamente, z es 
la posición relativa de las moléculas en la caja de simulación y z0 es una posición de 
referencia, normalmente la Superficie Divisoria de Gibbs (SDG) y d es un parámetro de 
ajuste atribuido al espesor de la interfase. Existen otras funciones para describir el perfil de 
densidades en equilibrios líquido–líquido, similares a la ecuación 2.29 para cada componente 
[31]. 
A A
A
B ABe
1
( ) rf
22 c
z z
z
t
  
 
  
 
 (2.30a) 
O O
O
B OBe
1
( ) rf
22 c
z z
z
t
  
 
  
 
 (2.30b) 
donde A y O son la densidad de la fase acuosa y orgánica, respectivamente, el subíndice B 
indica que es el valor promedio en el bulto de la respectiva fase, zA y zO son las posiciones 
promedio de la superficie de Gibbs de cada sustancia en la interfase y tc es un parámetro 
relacionado con el tamaño de la interfase en un equilibrio líquido–líquido, denominado 
criterio 10–50, propuesto por Senapati y Berkowitz [32]. Para determinar la posición de la 
Superficie Divisora de Gibbs (SDG) se usó el criterio equimolar, i.e., no hay acumulación 
del componente de referencia en la interfase. 
0
0
( ) ( )
z
z
i i i i id zz dz z
    


           (2.31) 
donde i es la adsorción del componente i en la superficie, z0 es la posición en la i=0 para 
el componente de referencia y los superíndices indican el valor en el bulto de las fases  y . 
Es común fijar la posición de la SDG como z=0. La Ec. 2.31 se evaluó numéricamente. 
2.2.2 Tensión superficial 
La tensión superficial e interfacial se determinó a partir de la relación 
32 
 
 
0
1
( ) ( )
2
2
zL
P N T
z
N T
P z P z dz
L
P P
  
    

 (2.32) 
donde Lz es la longitud de la caja de simulación, PN y PT son los componentes normal y 
tangencial del tensor presión, respectivamente. El tensor presión se determinó a partir del 
teorema del virial 
 cin
2
V
P E  (2.33) 
donde Ecin y  son la energía cinética y el virial en forma tensorial. En sistemas líquido–
vapor es necesario corregir la tensión superficial debido al truncamiento de las interacciones 
de no–enlace, por lo que se adiciona la de manera similar a la corrección en la presión. Esta 
corrección es ~3 mN m–1, mientras que para sistemas líquido–líquido este valor corresponde 
a ~0.5 mN m–1. El término de corrección se realizó siguiendo la metodología propuesta por 
Lundberg y Edholm [33], la cual asume que el perfil de densidad sigue una función 
hiperbólica. 
2.2.3 Parámetro de orden orientacional 
La orientación promedio de las moléculas se determinó vía el parámetro de orden 
orientacional, el cual está relacionado con el ángulo formado entre un vector característico 
de la molécula y el vector normal a la interfase. El parámetro de orden orientacional en 
términos de los polinomios de Legendre es 
0
2( ) 2
1
n
n k
n
k
n k
P
k
x x
n
k
 
  

 




 (2.34) 
donde n es el grado del polinomio. En este trabajo su valor es definido por la forma de la 
molécula y su vector característico. La molécula de agua, que presenta una estructura 
simétrica, puede ser descrita en término de vectores característicos direccionales, i.e., 
momento dipolar a, enlace OH b y vector normal al plano n⊥ (HOH) c, con grado n=1. 
Mientras que, para hidrocarburos lineales simétricos con extremos indistintos (caracterizados 
por un vector definido como el vector promedio que une a todos pseudo–átomos que 
conforman la molécula) n=2. 
1 cosP  (2.35) 
 
a 
i ii
q r 
b OH H O r rd 
c 
1 2OH OH
 n dd 
33 
 
2
2
1
3cos 1
2
zSP    (2.36) 
donde  es el ángulo formado entre el eje z y el vector característico. El valor del parámetro 
de orden orientacional puede variar según sea el valor de n. Para n=1 el parámetro de orden 
orientacionalpuede tomar valores entre 1 y –1 que corresponden a orientaciones 
perpendiculares a la interfase con direcciones opuestas, mientras que valores nulos sugieren 
una orientación preferencial paralela a la interfase o isotropía (Figura 2.7). 
 
 
Figura 2.7 Representación del parámetro de orden orientacional con vectores orientados. a) 
Vector perpendicular a superficie con P1=1 y P1= –1 o P2=1; b) Vector paralelo a superficie 
con P1=0 o P2= –0.5; c) Orientación de conjunto isotrópica con P1=0 o P2=0. 
En el caso de los alcanos, el parámetro de orden orientacional (P2 o Sz) puede variar entre –
0.5 y 1. Un valor promedio negativo sugiere que las moléculas se orientan preferentemente 
paralelas a la interfase, mientras que valores positivos sugieren que las orientaciones 
preferentes son perpendiculares a la interfase. Valores nulos denotan isotropía por parte de 
las moléculas. De manera similar al cálculo de densidades locales, el parámetro de orden 
orientacional se obtiene en segmentos con espesor de 0.5 nm. 
2.2.4 Distancia de extremo a extremo y radio de giro 
Las conformaciones moleculares se pueden inferir analizando la distancia de extremo a 
extremo de las moléculas (RF=||rC1–rCN||) y su radio de giro (Rg), los cuales están relacionados 
c) 
34 
 
al tamaño de moléculas con forma de cadena, e.g., polímeros, hidrocarburos lineales, etc. 
[34]. Así como su tendencia a elongarse o retraerse por su “elasticidad”. Se define el radio 
de giro como 
1
2
2
1
2 2 2
com
1
g i i
i
R m
M
  rr (2.37) 
donde M es la masa total de la molécula, mi es la masa del pseudo–átomo, ri y rcom son la 
posición del pseudo–átomo i y el centro de masa, respectivamente. Los componentes del 
radio de giro pueden proveer de una perspectiva de la simetría del sistema, los cuales se 
definen como 
1
2
1
2
2 2 2
com
1
g i ij j
j i
R m
M


 
 
 
   rr (2.38) 
donde  y j son las direcciones x, y o z. RF, Rg, Rgx, Rgy y Rgz son promedios determinados en 
segmentos con un espesor de 0.3 nm. 
2.2.5 Coeficiente de difusión 
El coeficiente de difusión es una propiedad de transporte macroscópica fácilmente medible 
vía DM. Sin embargo, para sistemas no–homogéneos, la ecuación de desplazamiento 
cuadrático medio (relación de Einstein) no es válida. Liu et al. [35] y Wick y Dang [36] 
presentan un método para calcular el coeficiente de difusión paralelo como función de la 
posición. Su argumento es que, si los segmentos que dividen la caja son paralelos a la 
interfase y lo suficientemente angostos, el coeficiente de difusión no variará respecto a z 
dentro del segmento. La difusión perpendicular se puede obtener resolviendo la ecuación de 
Smoluchowski o haciendo simulaciones de dinámica de Langevin. En este trabajo se siguió 
la metodología propuesta para el cálculo de coeficientes de difusión paralelos a la interfase 
(Dxy), determinados de la siguiente forma 
2( )
( ) lim
4 ( )
k
xy k
S
D


 


r
 (2.39) 
donde  es el intervalo de tiempo de la medición, S() es la probabilidad de permanencia de 
las moléculas en el segmento k durante el intervalo de medición y r2 es el Desplazamiento 
Cuadrático Medio (DCM) de la molécula i [37]. El DCM se calcula con la ecuación 
22 1( ) ( ) (0)
(0)
i ik
i kkN
 

  r r r (2.40) 
donde los paréntesis angulares con el subíndice k indican un promedio temporal en el 
segmento k, Nk(0) es el número de moléculas en el segmento k al inicio del intervalo de 
medición. Los coeficientes de difusión son resultado de la regresión lineal del tiempo sobre 
el DCM en la región lineal. Los coeficientes de difusión se calcularon en segmentos con 
espesor de 0.3 nm. 
35 
 
2.2.6 Energía libre 
Acorde a la TGIF, la tensión superficial y los perfiles de la energía libre están relacionados 
(véase Ec. 2.44). Las energías libres se calcularon usando un muestreo de sombrilla, como es 
descrito por Ghoufi et al. [38] y Biscay et al. [39]. En este trabajo se obtuvieron los perfiles 
de energía libre vía muestreo de sombrilla y el Método de Análisis de Histograma Ponderado, 
en inglés WHAM [40,41]. El cambio de energía libre F(z) se define como 
0
( )
( ) ln
( )
B
P z
F z k T
P z
   (2.41) 
donde F está en función de la posición de la molécula a lo largo del eje z y P(z) es la función 
de distribución promedio 
 
( )
exp
( )
( )
w
w
N
ii
N
j j ji
h z
n w
P
z
z
f 
 




 (2.42a) 
 exp exp ( ) ( )j jf w z P z dz      (2.42b) 
donde hi es el histograma generado, Nw es el número total de puntos generados en el 
histograma,  es el inverso de la temperatura 1/kBT, wj es el potencial armónico entre la 
molécula y una posición de referencia y fj es una constante de energía libre. Para cada 
posición j (ventana) se restringe una molécula sonda al potencial 
 
2
( )
2
ci
i iw
K
z z z  (2.43) 
donde zic es una posición a la que se restringe la molécula y Ki es la constante de fuerza, la 
cual se fijó al valor 1000 kJ mol–1nm–2. Para generar las configuraciones iniciales se llevó a 
cabo una simulación guiada, tomando como referencia a una molécula (A) ubicada a la mitad 
de la caja de simulación y una molécula sonda (B) a una distancia fija de la molécula A. La 
molécula B se aleja de la molécula A a lo largo del eje z desde el bulto de la fase líquida hacia 
el bulto de la fase vapor con incrementos de 0.2 nm. Las ventanas generadas se usan como 
configuraciones iniciales para el muestreo de simulaciones de sombrilla. 
 
 
2.3 Teoría gradiente 
 
La Teoría Gradiente de Interfases Fluidas (TGIF) expresa la densidad de la energía de 
Helmholtz de un fluido no–homogéneo como una serie de Taylor truncada al segundo 
término centrada en la energía de Helmholtz de un fluido homogéneo evaluada a la densidad 
local. En la TGIF se asume que los gradientes moleculares en la interfase son pequeños 
comparados con el valor reciproco de la distancia intermolecular, de manera que la densidad 
y sus derivadas se pueden tratar como variables independientes. También considera que la 
36 
 
energía de Helmholtz tiene dos contribuciones, una que corresponde al fluido homogéneo y 
otra debido a las no–homogeneidades en la interfase. 
 
2
0
1
( )
2
d
F A f z c dz
dz




  
   
   
 (2.44) 
donde A es el área interfacial, f0 es la densidad de energía de Helmholtz del fluido 
homogéneo,  es la densidad local, z es la coordenada normal a la interfase y c es el parámetro 
de influencia, el cual contiene información acerca de la estructura molecular y la respuesta 
del gradiente de la densidad a cambios del potencial químico local de su valor promedio en 
el bulto de las fases fluidas. Usualmente, se considera a c como un parámetro que no depende 
de la densidad local, de manera que la ecuación de Euler–Lagrange que minimiza la energía 
de Helmholtz es 
2
2
d
c
dz







 (2.45) 
donde  es el gran potencial y =f()–+P. De manera que el perfil de densidad, 
solución de la ecuación 2.45, así como la tensión superficial se calculan por la siguiente 
ecuación 
2
d
c dz
dz




 
  
 
 (2.46) 
Las simulaciones moleculares pueden generar perfiles de densidad y tensiones superficiales, 
de manera que la única variable desconocida en la ecuación 2.46 es el parámetro de 
influencia. Por otro lado, la combinación de las ecuaciones 2.45 y 2.46 permite calcular la 
tensión superficial vía Ecuaciones de Estado (EdE): 
 
1
22c d




    (2.47) 
La energía de Helmholtz del fluido homogéneo y sus potenciales químicos se determinan 
usando cualquier EdE. En este trabajo se usaron las ecuaciones SRK y PC–SAFT debido al 
éxito reportado en el cálculo de propiedades termodinámicas volumétricas, así como el 
cálculo de tensiones superficiales e interfaciales.Cabe mencionar que para la solución de la Ec. 2.45, así como el cálculo del gran potencial 
se requiere de las condiciones a la frontera, i.e., las propiedades de las fases fluidas en 
contacto en el equilibrio, i.e., equilibrio térmico T=T, equilibrio mecánico P=P y 
equilibrio material 1=1, …, N=N. 
 
2.4 Ecuaciones de estado 
 
37 
 
Las Ecuaciones de Estado son relaciones matemáticas que describen la dependencia de las 
variables termodinámicas intensivas 
 , ,f f T P v (2.48) 
Para el uso de la teoría gradiente se requiere conocer el comportamiento volumétrico al 
equilibrio, potenciales químicos y la energía libre de Helmholtz. Estas propiedades se 
determinaron usando dos EdE populares; PC–SAFT y SRK. 
2.4.1 Teoría estadística de fluidos asociantes (SAFT) 
El cálculo de propiedades de los denominados fluidos simples resulta una tarea que no 
presenta mayor problema debido a que son fluidos cuyas moléculas son casi esféricas y 
químicamente inertes. Sin embargo, fluidos que presentan una cadena larga o interacciones 
de específicas altamente direccionales y de corto alcance representan un gran problema. 
Uno de los modelos más exitosos capaz de describir este tipo de fluidos y sus interacciones 
es el desarrollado por Chapman et al. [42–45], quienes usando la teoría termodinámica de 
perturbación de primer orden de Wertheim dieron origen a la denominada Teoría Estadística 
de Fluidos Asociantes, SAFT. Wertheim obtuvo una expresión simplificada de la energía de 
Helmholtz y propuso una relación entre la energía de Helmholtz residual y la densidad del 
monómero debido al fenómeno de asociación. Dentro de las aproximaciones que involucran 
los modelos de la familia SAFT se encuentran: 
• Las estructuras ramificadas o cadenas están permitidas dentro de la teoría, siendo lo 
contrario para estructuras tipo anillo. 
• Solo se permiten enlaces sencillos por sitio de asociación. Las restricciones que 
considera esta aproximación son: dos núcleos de moléculas asociadas evitan que una 
tercera molécula se asocie por efectos de estericidad; no se puede asociar un sitio de una 
molécula con dos sitios de otra molécula; no se permiten enlaces dobles entre dos 
moléculas. 
• El ángulo formado entre los sitios de enlace es independiente de las propiedades del 
fluido. 
• La actividad de un sitio es independiente de la actividad de otros sitios de la misma 
molécula. 
Chapman et al. [42–45] y posteriormente Huang y Radosz [46,47] presentaron la EdE SAFT 
que considera a un fluido de referencia que considera tanto la forma e interacciones de las 
moléculas. El modelo SAFT considera la energía de Helmholtz residual como la suma de las 
contribuciones debidas a la interacción entre segmentos, la formación de cadenas del mismo 
tipo de segmentos y la formación de enlaces de asociación entre dos segmentos o cadenas: 
seg cd asocRF F F F   (2.49) 
donde FR, Fseg, Fcd y Fasoc son las energías de Helmholtz residual, de segmento, de cadena 
dura y de asociación, respectivamente (Figura 2.8). 
38 
 
Los monómeros o segmentos, son la unidad básica para la construcción de la molécula. Estos 
segmentos tienen la capacidad de interactuar a través de fuerzas de dispersión y repulsión. 
Esta interacción no–direccional da lugar a la formación de la molécula tipo cadena formada 
por m segmentos conectados unos a otros. La contribución a la energía de Helmholtz residual 
que corresponde a la contribución individual de todos los segmentos que forman la molécula 
está dada por 
segseg
0
i
B i
i
FF
m
Nk T N T
x
R
  (2.50) 
 
Figura 2.8 Representación de etapas de formación de cadenas y complejos. 
Donde F0
seg es la energía de Helmholtz residual de segmentos esféricos no asociados. Si Los 
segmentos se componen de esferas duras con partes dispersivas 
seg ed disp
0 0 0F F F  (2.51) 
donde el término de esferas duras se reemplaza por la EdE propuesta por Carnahan y Starling 
[48] 
2ed
0
2
4
(1 )
3
B
F
Nk T
 


 (2.52) 
3
6
i i i
i
m dx

   (2.53) 
donde  es la fracción de empaquetamiento o densidad reducida. Para mezclas las 
expresiones anteriores suponen que la energía de Helmholtz es la de una esfera dura junto 
con la teoría de un fluido propuesta por van der Waals para el término dispersivo [47]. 
Fluido de segmentos rígidos 
Interacción por 
asociación 
Interacción 
por atracción 
Formación de cadenas 
39 
 
Para la contribución de cadenas duras se considera que los segmentos se polimerizan al 
establecer contacto en el borde de la molécula, en el límite de fuerza de enlace infinito. 
e
cd
d(1 ) lni i
i
ij
B
F
x
N
gm
k T
  (2.54) 
donde gij
ed la función de distribución radial propuesta por Mansoori et al. [49] que depende 
del diámetro efectivo de colisión d(T). 
2
2
2 2
2 3
3 3
e
3
d 3 21
1 (1 ) (1 )
i j i j
i j i j
ij
d d
g
d
d dd d
d  
  
 
 
   
          
 (2.55) 
Cabe mencionar que la teoría TTP1 de Wertheim produce una buena aproximación para 
moléculas lineales; sin embargo, no considera de manera estricta los ángulos de enlace, lo 
cual puede hacerse al considerar ordenes de perturbación mayores, lo cual conlleva a 
expresiones complicadas e involucran funciones de correlación de tres o más cuerpos para el 
fluido de referencia. 
El término de asociación, que describe interacciónes de corto alcance altamente direccionales 
como puentes de hidrógeno o transferencia de cargas, es el corazón de las ecuaciones de la 
familia SAFT, siendo este el responsable de cuantificar interacciones como puentes de 
hidrógeno. A partir de la teoría de perturbación de Wertheim, Huang y Radosz propusieron 
el siguiente modelo de contribución a la energía de Helmholtz residual para fluidos que tienen 
la capacidad de formar puentes de hidrógeno (asociarse): 
asoc
ln
2 2
i
i
ii A
A
A i
i
B
MF
x
N T
X
k
X  
  

 
  
  (2.56) 
Donde xi es la fracción molar de la especie i, XAi es la fracción molar de moléculas i no 
asociadas al sitio A y Mi es el número de sitios de asociación por molécula. La fracción de 
moléculas no asociadas al sitio A es 
1
Av1
i j
j
ji
B A BA
j
j B
NX x X

 
   
 
 
 (2.57) 
donde NAv es el número de Avogadro,  es la densidad molar de las moléculas y AiBj es la 
fuerza de asociación, propiedad clave que caracteriza el enlace de asociación. 
3 ( ) exp 1
i j
i j i j
A B
A B A B
ij ij
k
g d
T
 
 
  



 

 (2.58) 
dondeAiBj es la energía de asociación y AiBj el volumen de asociación. La forma y el valor 
de las fuerzas de asociación (Ecs. 2.57 y 2.58) depende del esquema de asociación propuesto 
para la molécula. Dicho esquema se esquematiza por Huang y Radoz [46] en la Tabla VII. 
40 
 
Existen una gran variedad de ecuación pertenecientes a la familia SAFT. Entre ellas destaca 
la modificación realizada al término dispersivo por Gross y Sadowski [50], quienes 
extendieron la teoría de perturbación de Barker y Henderson [51,52] y tomaron como fluido 
de referencia a una cadena dura, a diferencia del modelo original, donde el fluido de 
referencia son esferas duras. 
cd disp asocRF F F F   (2.59) 
En el modelo PC–SAFT, la expresión para la energía de Helmholtz del término dispersivo es 
1
disp
2
B B B
A AF
Nk T Nk T Nk T
  (2.60) 
donde A1 representa al promedio de pares de segmentos en el rango de interacciones atractivo, 
mientras que A2 está relacionado con el momento de la distribución del promedio. 
 
2.4.2 Ecuación de estado SRK 
La ecuación cúbica de Soave–Redlich–Kwong, SRK, es una modificación propuesta por 
Soave [53] a la EdE Redlich–Kwong. Esta ecuación considera la dependencia con la 
temperatura del término atractivo. Soave propuso una dependencia con el factor de acéntrico