Berechnung von Grenzflächeneigenschaften mithilfe der klassischen Dichtefunktionaltheorie mit einem Funktional auf Grundlage der PCP-SAFT Zustandsgleichung Von der Fakultät 4 – Energie-, Verfahrens- und Biotechnik – der Universität Stuttgart zur Erlangung der Würde eines Doktors der Ingenieurwissenschaften (Dr.-Ing.) genehmigte Abhandlung Vorgelegt von M.Sc. Christoph Klink aus Herrenberg Hauptberichter: Prof. Dr.-Ing. J. Gross Mitberichter: Prof. Dr. habil. rer. nat. S. Enders Tag der mündlichen Prüfung: 19.11.2018 Institut für Technische Thermodynamik und Thermische Verfahrenstechnik der Universität Stuttgart 2018 III Inhaltsverzeichnis Symbolverzeichnis V Kurzfassung XI Abstract XIII 1 Einleitung 1 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Zielsetzung der Arbeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 2 Thermodynamische Grundlagen 4 2.1 Beschreibung der Phasengrenzfläche . . . . . . . . . . . . . . . . . . . . . 4 2.2 Die Helmholtzenergie und ihre Ableitungen . . . . . . . . . . . . . . . . . 6 3 PCP-SAFT Zustandsgleichung für Mischungen 9 3.1 Molekulares Modell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 3.2 Wechselwirkungsbeiträge . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 3.2.1 Hart-Kugel-Beitrag und Beitrag harter Ketten . . . . . . . . . . . 10 3.2.2 Beitrag durch dispersive Wechselwirkungen . . . . . . . . . . . . . 11 3.2.3 Beitrag durch assoziierende Wechselwirkungen . . . . . . . . . . . 12 3.2.4 Beitrag durch multipolare Wechselwirkungen . . . . . . . . . . . . 13 4 Dichtefunktionaltheorie 16 4.1 Das Großkanonische Potential Ω . . . . . . . . . . . . . . . . . . . . . . . 16 4.2 Funktionalbeiträge zur Helmholtzenergie . . . . . . . . . . . . . . . . . . 17 4.2.1 Idealer Gasbeitrag . . . . . . . . . . . . . . . . . . . . . . . . . . 18 4.2.2 Hart-Kugel-Beitrag - Fundamental Measure Theorie . . . . . . . . 19 4.2.3 Kettenbeitrag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.2.4 Entwicklung eines neuen Helmholtzenergiefunktionals dispersiver Wechselwirkungen für Mischungen . . . . . . . . . . . . . . . . . 21 4.2.5 Entwicklung eines neuen Helmholtzenergiefunktionals dipolarer Wechselwirkungen . . . . . . . . . . . . . . . . . . . . . . . . . . 24 IV 5 Grundlagen der Nichtgleichgewichts-Thermodynamik an Grenzflächen 29 5.1 Entropieproduktion aus Wärme- und Stofftransport durch Grenzflächen . 29 5.2 Integralbeziehungen zur Beschreibung von Transportwiderstandsko- effizienten an Grenzflächen . . . . . . . . . . . . . . . . . . . . . . . . . . 32 6 Berechnungsablauf und numerische Umsetzung 36 6.1 Grenzflächenspannung dipolarer Reinstoffsysteme . . . . . . . . . . . . . 36 6.2 Grenzflächenspannung von binären Dampf-Flüssig-Systemen . . . . . . . 40 6.3 Grenzflächenspannung von binären Flüssig-Flüssig-Systemen . . . . . . . 41 6.4 Transportwiderstände . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.4.1 Reinstoffsysteme . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 6.4.2 Multikomponentenmischungen . . . . . . . . . . . . . . . . . . . . 44 7 Ergebnisse 47 7.1 Grenzflächeneigenschaften dipolarer Reinstoffe . . . . . . . . . . . . . . . 47 7.2 Grenzflächeneigenschaften binärer Dampf-Flüssig-Systeme . . . . . . . . 52 7.3 Grenzflächeneigenschaften binärer Flüssig-Flüssig-Systeme . . . . . . . . 59 7.4 Transportwiderstände . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 7.4.1 Reinstoffe . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 7.4.2 Binäre Mischung . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 8 Zusammenfassung 92 9 Ausblick 94 Anhang 96 A Analyse nicht-winkelabhängiger Helmholtzenergiebeiträge 97 B Partielle molare Helmholtzenergiebeiträge 98 B.1 Hart-Kugel-Beitrag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 B.2 Kettenbeitrag . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99 B.3 Beitrag durch dispersive Wechselwirkungen . . . . . . . . . . . . . . . . . 100 Literaturverzeichnis 101 V Symbolverzeichnis Lateinische Symbole Symbol Einheit Bedeutung A [m2] Fläche ai [-] Modellkonstanten im Dispersionsterm der PCP- SAFT Zustandsgleichung an,ij [-] Modellkonstanten im multipolaren Wechselwir- kungsterm der PCP-SAFT Zustandsgleichung bi [-] Modellkonstanten im Dispersionsterm der PCP- SAFT Zustandsgleichung bn,ij [-] Modellkonstanten im multipolaren Wechselwir- kungsterm der PCP-SAFT Zustandsgleichung cn,ij [-] Modellkonstanten im multipolaren Wechselwir- kungsterm der PCP-SAFT Zustandsgleichung D [m2/s] Diffusionskoeffizient d(T ) [Å] Temperaturabhängiger Segmentdurchmesser F [J] Helmholtzenergie; Freie Energie f [bar] Fugazität [J/m3] Helmholtzenergiedichte g(r) [-] Radiale Paarverteilungsfunktion H [J] Enthalpie h; hi [J/mol] Molare bzw. partielle molare Enthalpie i; j; k [-] Laufvariablen I1; I2 [-] Parameter im Dispersionsterm der PCP-SAFT Zu- standsgleichung J [mol/s] Strom J2; J3 [-] Parameter im dipolaren Wechselwirkungsterm der PCP-SAFT Zustandsgleichung VI Symbol Einheit Bedeutung k [J/K] Boltzmann-Konstante k = 1, 380 · 10−23 J/K kij [-] Binärparameter lij [-] Binärparameter m [-] Segmentanzahl m(ω) [-] Einheitsvektor der Orientierungswinkel des Dipol- Dipol-Potentials N [mol] Gesamtmolekülzahl Ni [mol] Molekülzahl einer Komponente i n [mol] Gesamtstoffmenge [-] Anzahl an Komponenten in einer Mischung ni [mol] Stoffmenge einer Komponente i NAv [mol−1] Avogadro-Konstante NAv = 6, 022 · 1023 mol−1 p [bar] Druck Q [D · Å] Quadrupolmoment R [J/(mol K)] Universelle Gaskonstante R = 8, 314 J/mol ·K Rqq [(m2 s)/(J K)] Widerstandskoeffizient gegen reinen Wärmetrans- port Rqi [(m2 s)/(mol K)] Widerstandskoeffizient gegen gekoppelten Wärme- und Stofftransport Rii [(J m2 s)/(mol2 K)] Widerstandskoeffizient gegen reinen Stofftransport r [m] Abstands- bzw. Ortsvariable r(z) [(m s)/(J K)] Lokaler Widerstandskoeffizient gegen reinen Wär- metransport [(m s)/(mol K)] Lokaler Widerstandskoeffizient gegen gekoppelten Wärme- und Stofftransport [(J m s)/(mol2 K)] Lokaler Widerstandskoeffizient gegen reinen Stoff- transport S [J/K] Entropie s; si [J/(mol K)] Molare bzw. partielle molare Entropie ST [K−1] Soretkoeffizient T [K] Temperatur VII Symbol Einheit Bedeutung uDD [J] Dipol-Dipol-Potential V [m3] Volumen v; vi [m3/mol] Molares bzw. partielles molares Volumen X [J/(mol2 K m s)] Konjugierte Triebkraft xi [mol/mol] Stoffmengenanteil XAi [-] Stoffmengenanteil einer Komponente i im assoziie- renden Wechselwirkungsanteil der PCP-SAFT Zu- standsgleichung, der nicht an einer Bindungsstelle A gebunden ist XBi [-] Stoffmengenanteil einer Komponente i im assoziie- renden Wechselwirkungsanteil der PCP-SAFT Zu- standsgleichung, der nicht an einer Bindungsstelle B gebunden ist yddii [-] Kavitätskorrelationsfunktion Z [-] Allgemeine Zustandsgröße [-] Kompressibilitätsfaktor z [m] Abstands- bzw. Ortsvariable zi [-] allgemeine partielle molare Zustandsgröße Griechische Symbole Symbol Einheit Bedeutung α [(m9 s)/(mol2 J K)] anpassbarer Parameter im lokalen Widerstands- profil gegen Wärmetransport [(m9 s)/(mol3 K)] anpassbarer Parameter im lokalen Widerstands- profil gegen gekoppelten Wärme- und Stofftrans- port α(z) [-] Orientierungsverteilungsfunktion β [-] anpassbarer Parameter im lokalen Widerstands- profil [J−1] reziproke Temperatur β = 1 kT γ [N/m] Grenzflächenspannung VIII Symbol Einheit Bedeutung ∆AiBi [-] Assoziierungsstärke δ [-] Dirac-Funktion ε [J] Energieparameter εAiBi [J] Assoziierungsenergie η [-] Packungsdichte Θ [-] Heaviside Sprungfunktion κ [-] Kriterium für die Grenzflächendicke [-] Kapillarwellenverhältnis κAiBi [-] Assoziierungsvolumen Λ [Å] DeBroglie-Wellenlänge λ [W/(m K)] Wärmeleitfähigkeitskoeffizient λi [-] Gemittelte Dichte am Berührungspunkt um ein Segment der Kette i µ [J/mol] Chemisches Potential π [-] Koexistierende Phase ρ [Å−3] Teilchendichte ρ̂ [Å−3] Erweiterte Teilchendichte σ [Å] Segmentdurchmesser [J/(mol K)] Entropieproduktion σ(z) [J/(mol K m)] Lokale Entropieproduktion ϕ [-] Fugazitätskeoffizient φ [J] Potential χ [-] Akzeptanzkriterium der Dichteiteration Ω [J] Großkanonisches Potential ω [-] Ausrichtungsvektor des Dipolmoments IX Indizes (hochgestellt) ′, ′′, π Phase θ Referenzpunkt 1PT Störungstheorie erster Ordnung assoc Assoziierend dd, EoS Dipol-Dipol bezogen auf die Zustandsgleichung DD Dipol-Dipol disp Dispersiv e Bezogen auf die Gesamtenergie g Bezogen auf die Gasphase hc Harte-Kette (engl. hard chain) hs Harte-Kugel (engl. hard sphere) ig; id. Gas Ideales Gas iso Isotrop l Bezogen auf die Flüssigphase L1; L2 Flüssigphase L1 bzw. L2 PT Störungstheorie QQ Quadrupol-Quadrupol rein Reinstoff res Residuell S Siede- Indizes (tiefgestellt) bulk Bezogen auf eine Kernphase c Bezogen auf den kritischen Punkt DFT Dichtefunktionaltheorie e Bezogen auf den Gesamtenergiefluss q Bezogen auf den Wärmefluss qi Bezogen auf den gekoppelten Wärme- und Stofffluss ξi Bezogen auf den Stofffluss X Abkürzungen BMCSL Booblik-Mansoori-Carnahan-Starling-Leland DDK Dampfdruckkurve DFT Dichtefunktionaltheorie DMF Dimethylformamid FMT Fundamental measure theory G Gasphase KP Kritischer Punkt L Flüssigphase LLE Flüssig-Flüssig-Gleichgewicht (engl. liquid liquid equilibrium) MAE Mittlerer absoluter Fehler (engl. mean absolute error) MAPE Mittlerer absoluter prozentualer Fehler (engl. mean absolute percentage error) MeOH Methanol NEMD Nichtgleichgewichtsmolekulardynamik (engl. non-equilibrium molecular dynamics PCP − SAFT Perturbed chain polar statistical associating fluid theory S Feststoffphase SAFT − γ Als Gruppenbeitragsmethode formulierte statistical associa- ting fluid theory SAFT − V R Statistical associating fluid theory with variable range Schmelz −DK Schmelzdruckkurve SGT Dichtegradiententheori (engl. square gradient theory) Sub.−DK Sublimationsdruckkurve THF Tetrahydrofuran TP Tripelpunkt UCST Obere kritische Entmischungstemperatur (engl. upper critical solution temperature) V LE Dampf-Flüssig-Gleichgewicht (engl. vapor liquid equilibrium) XI Kurzfassung Diese Arbeit untersucht die Vorhersage von Grenzflächeneigenschaften fluider Mischun- gen. Um Größen wie die Grenzflächenspannung oder Widerstandskoeffizienten einer Grenzfläche gegenüber Wärme- und Stofftransport berechnen zu können, ist – neben einer möglichst genauen Zustandsgleichung zur Bestimmung der Kernphasengrößen – ein Formalismus zur Berechnung des Dichteverlaufs über die Grenzfläche hinweg erfor- derlich. Ein Funktional der Helmholtzenergie auf Basis der PCP-SAFT Zustandsgleichung ist in der Literatur bereits gegeben. Dieses wurde – eingebettet in einen Dichtefunktionaltheo- rie (DFT) Formalismus – zur Berechnung der Grenzflächenspannung von Dampf-Flüssig- Gleichgewichten nicht oder schwach polarer Reinstoffe mit sehr hoher Genauigkeit im Vergleich zu experimentellen Daten verwendet. Ziel dieser Arbeit war zunächst, den vorgestellten DFT Formalismus für stark dipolare Reinstoffe und schwach dipolare Mi- schungen zu erweitern. Für die Betrachtung dipolarer Reinstoffe war es erforderlich, den gesamten Formalismus winkelabhängig zu formulieren, um die Ausrichtung des Dipol- moments berücksichtigen zu können. Der zusätzliche Wechselwirkkungsbeitrag im Helm- holtzenergiefunktional des DFT Formalismus wurde durch eine Störungstheorie erster Ordnung beschrieben. Der erwünschte Effekt, die Grenzflächenspannung im Vergleich zum nicht-polaren DFT Formalismus durch Freigabe des zusätzlichen Freiheitsgrads der Orientierung zu erniedrigen, konnte nur in schwachem Maße erreicht werden. Für die Berechnung von Mischungen wurde das Helmholtzenergiefunktional für die di- spersiven Wechselwirkungen auf Basis einer Drei-Fluid-Theorie entwickelt. Die Grenz- flächenspannung der Dampf-Flüssig-Grenzfläche einer Reihe nicht- oder schwach pola- rer, binärer Mischungen konnte mit hoher Genauigkeit im Vergleich zu experimentellen Daten berechnet werden. Hierzu war lediglich ein Parameter für die Beschreibung der Mischungseffekte zwischen strukturell unterschiedlichen Komponenten notwendig, der an experimentelle Daten des Dampf-Flüssig-Gleichgewichts angepasst wurde. Es wur- de darauf verzichtet, Parameter an explizite Zielgrößen wie die Grenzflächenspannung anzupassen, um den prädiktiven Charakter des Formalismus beizubehalten. Die Be- rechnung von binären Flüssig-Flüssig-Gleichgewichten zeigte ebenfalls unerwartet gute Ergebnisse für die Grenzflächenspannung, sofern das Phasenverhalten durch die verwen- dete Zustandsgleichung gut beschrieben wurde. XII Um die Vielseitigkeit des weiterentwickelten Formalismus in einem weiteren Kontext zu prüfen, wurde er angewandt, um die Widerstandskoeffizienten der Grenzfläche gegen- über gekoppeltem Wärme- und Stofftransport zu berechnen. Diese Größen waren bisher lediglich mittels Molekularsimulationen zugänglich, welche nur unter hohem Rechenauf- wand durchführbar sind. Für Reinstoffe ließ sich der Formalismus mit sehr guter Über- einstimmung mit Daten aus Molekularsimualtionen anwenden. Aufgrund der schlechten Datenbasis der Widerstandskoeffizienten für Mischungen konnte die Theorie nur für eine Modellmischung überprüft werden, wobei sich ebenfalls eine sehr gute Übereinstimmun- gen zeigte. Auf Basis dieser Rechnungen ist die Erwartung gegeben, dass auch für reale Mischungen gute Ergebnisse erzielt werden. XIII Abstract This study investigates the prediction of interfacial properties of fluid mixtures. To be able to calculate quantities like the interfacial tension or resistivities of an interface against thermal or mass transport, a formalism – besides an equation of state as precise as possible to determine bulk quantities – is needed to calculate the density development across the interface. A Helmholtz energy functional based on the PCP-SAFT equation of state is already known in literature. This was applied – embedded in a density functional theory (DFT) formalism – to calculate the interfacial tension of vapor-liquid equilibria of pure non- or weakly polar substances with very high accuracy compared to experimental data. Initi- ally, the aim of this study was to extend the present DFT formalism to strong dipolar substances and weakly dipolar mixtures. The examination of dipolar substances made it necessary to formulate the complete formalism in an angle dependent way to consi- der the orientation of the dipole moment. The additional contribution to the Helmholtz energy due to angle dependent interactions was described by a perturbation theorie of first order. The desired effect to lower the interfacial tension by releasing the orientation as additional degree of freedom compared to the non-polar DFT formalism could only be realized on a low degree of success. For the calculation of mixtures, the Helmholtz energy functional for dispersive interac- tions was developed based on a three-fluid-theory. The interfacial tension of the vapor- liquid interface of a series of non- or weakly polar binary mixtures could be calculated with high accuracy compared to experimental data. For this, only one parameter was necessary to describe the mixing effects between structurally different components which was adjusted to experimaental data of the vapor-liquid equilibrium. To retain the predic- tive character of the formalism, no parameter was adjusted to target quantities like the interfacial tension. The calculation of binary liquid-liquid euqilibria also showed unex- pectedly good results for the interfacial tension in case the phase behavior was described properly by the applied equation of state. To investigate the versatility of the refined formalism in a different context, it was ap- plied to calculate resistivities of the interface against coupled heat and mass transport. These quantities could be accessed so far only via molecular simulations which require a high amount of computing effort. For pure substances the formalism could be applied XIV with data out of molecular simulations with very high accuracy. Due to a poor data basis for transport resistivities for mixtures the theory could only be tested for a model mixture at which very good agreement resulted as well. Based on these calculations good results even for real mixtures are expected. 1 1 Einleitung 1.1 Motivation Thermodynamische Prozesse, bei denen Phasenübergänge in Mischungen stattfinden, erfordern eine genaue Kenntnis der Stoffdaten sowohl der beteiligten Reinstoffe als auch der Mischungseigenschaften. Neben einer Vielzahl an empirisch entwickelten Korrelatio- nen liefern Zustandsgleichungen auf Basis von Theorien aus der statistischen Mechanik akkurate Ergebnisse für den Zusammenhang von Druck, Volumen, Zusammensetzung und Temperatur. Mitte der 80er Jahre legte Wertheim1,2,3,4 mit seiner Graphentheorie zu kurzreichweitig anziehenden Wechselwirkungen den Grundstein für eine Familie von Zustandsgleichungen, den sogenannten SAFT-Zustandsgleichungen. Die statitical asso- ciating fluid Theorie (SAFT) abstrahiert Moleküle als Ketten aus gebundenen sphä- rischen Kugelsegmenten, die durch verschiedene Wechselwirkungsbeiträge miteinander interagieren5,6,7. Zur Familie der SAFT-Modelle zählen soft-SAFT8,9, SAFT-VR10,11, SAFT-γ 12,13 sowie PC-SAFT14,15 und die polare Version PCP-SAFT16,17,18,19. Detaillier- te Beschreibungen, Vergleiche und Anwendungsgebiete sind in den Arbeiten von Müller and Gubbins20, Economou21, Paricaud et al.22 und Tan et al.23 zusammengefasst. Alle SAFT Zusatandsgleichungen sind als Funktion von Temperatur und Dichte aller be- trachteten Komponenten in der Helmholtzenergie formuliert. Für Prozesse wie die Koaleszenz von Tropfen, bei denen zusätzlich zum Zusammenhang von Druck, Temperatur und Volumen auch Größen wie die Grenzflächenspannung ent- scheidend sind, werden weitere Formalismen und Berechnungsmodelle notwendig. Eine zentrale Eigenschaft einer Grenzfläche im Gleichgewicht ist ein kontinuierlicher Ver- lauf der Dichte (auf molekularer Längenskala) jeder betrachteten Komponente einer Mischung als Funktion des Ortes, welcher die Grundlage für die Berechnung weiterer Grenzflächeneigenschaften bildet. Neben Molekularsimulationen, die einen hohen Re- chenaufwand erfordern, haben sich zwei Theorien zur Berechnung des Dichteverlaufs als vielversprechend gezeigt. Die Dichte-Gradienten-Theorie (engl. square gradient theory, SGT) wurde von van der Waals und Rayleigh abgeleitet und von Evans24 maßgeblich weiterentwickelt. Hierbei wird die Helmholtzenergie als Reihenentwicklung des Dich- teprofils um den homogenen Verlauf entwickelt, die nach dem Term zweiter Ordnung – dem quadrierten Gradienten – abgeschnitten wird. Bei der Formulierung kommt in 2 KAPITEL 1. EINLEITUNG der Regel ein sogennanter Einflussparameter zum Einsatz, der empirisch an Reinstoff- oder Mischungseigenschaften angepasst wird. Die SGT wurde in zahlreichen Arbeiten mit verschiedenen Zustandsgleichungen kombiniert. Neben der Anwendung von kubi- schen Zustandsgleichungen wurden in diversen Arbeiten auch Varianten von SAFT- Zustandsgleichungen genutzt25,26. Die zweite Theorie zur Berechnung des Dichteverlaufs ist die Dichtefunktionaltheorie (DFT). In dieser Theorie wird das großkanonische Potential in allgemeiner Form als Funktional der Dichte zu gegebenem chemischen Potential µ, Volumen V und Tempe- ratur T formuliert. Verschiedene Vereinfachungen in der Formulierung des Helmholtz- energiefunktionals kommen zur Anwendung. In den meisten DFT Arbeiten wird eine Störungstheorie angewandt, bei der das Wechselwirkungspotential in einen repulsiven und einen attraktiven Anteil aufgespalten wird. Der attraktive Anteil wird hierbei als Störung betrachtet. Diese Unterteilung stellt an sich noch keine Näherung dar. Vielmehr die Wahl des Modells für das Referenzfluid sowie die Entwicklungen der abgeschnittenen Anteile der attraktiven Wechselwirkungen nutzen Vereinfachungen. Für Reinstoffe zeig- te Gross27 die erfolgreiche Anwendung der DFT in Kombination mit der PCP-SAFT Zustandsgleichung für nicht- oder schwach polare Reinstoffe. 1.2 Zielsetzung der Arbeit Das Ziel dieser Arbeit ist die Erweiterung des von Gross27 entwickelten DFT Formalis- mus sowohl auf stark polare Reinstoffe als auch auf Mischungen unterschiedlicher, nicht oder schwach dipolarer Komponenten. Als Zustandsgleichung kommt die PCP-SAFT Zustandsgleichung zum Einsatz. Für die Beschreibung stark polarer Reinstoffe wie beispielsweise Aceton ist eine Erwei- terung des bestehenden Formalismus unter Berücksichtigung der Winkelabhängigkeit des vorherrschenden Dipolmoments notwendig. Damit wird dem System ein zusätzli- cher Freiheitsgrad bereitgestellt, um seine freie Energie im Gleichgewicht zu minimieren. Außerdem muss ein neues Helmholtzenergiefunktional zur Beschreibung der dipolaren Wechselwirkungen entwickelt werden. Das Ziel ist die genauere Berechnung der Grenz- flächenspannung im Vergleich zu den Ergebnissen von Gross27. Die Beschreibung von Mehrkomponentensystemen erfordert ein neues Helmholtzener- giefunktional für die dispersiven Wechselwirkungen. Hier soll, wie durch Groß bereits für Reinstoffe formuliert, das neue Funktional die radiale Paarverteilungsfunktion be- rücksichtigen und weiterhin mit der PCP-SAFT Zustandsgleichung kompatibel sein. Es sollen sowohl Dampf-Flüssig-Gleichgewichte als auch Flüssig-Flüssig-Gleichgewichte binärer Gemische betrachtet werden und deren Grenzflächenspannungen als Zielgröße berechnet werden. KAPITEL 1. EINLEITUNG 3 Um den Anwendungsbereich des weiterentwickelten Formalismus für binäre Gemische über statische Grenzflächeneigenschaften wie die Grenzflächenspannung hinaus zu über- prüfen, werden dynamische Stoff- und Wärmetransportprozesse betrachtet. Hierbei steht die Berechnung von gekoppelten Wärme- und Stofftransportwiderständen binärer Dampf-Flüssig-Systeme über die Grenzfläche hinweg als weiteres Ziel im Fokus. Alle berechneten Größen werden mit experimentellen Daten oder Daten aus Molekular- simulationen aus der Literatur verglichen, um die Vorhersagen des Modells zu bewerten. 4 2 Thermodynamische Grundlagen 2.1 Beschreibung der Phasengrenzfläche Eine Grenzfläche (oder auch Phasengrenze) beschreibt im Allgemeinen die Fläche zwi- schen zwei benachbarten Phasen. Sowohl für Reinstoffsysteme als auch für Multikom- ponentenmischungen existiert eine Grenzfläche zwischen Phasen ungleichen Aggregat- zustandes wie fest und flüssig, flüssig und gasförmig oder fest und gasförmig. Zusätzlich können Grenzflächen zwischen Phasen gleichen Aggregatzustandes auftreten. Abbildung 2.1 zeigt schematisch den Dichteverlauf an einer Dampf-Flüssig-Grenzfläche. Im Gegensatz zur kontinuumsmechanischen Annahme einer sprunghaften Änderung der Dichte an der Grenzfläche (gestrichelte Linie), ergibt sich ein kontinuierlicher Verlauf zwischen den Kernphasenwerten (durchgezogene Linie). Im thermodynamischen Gleichgewicht gelten drei Gleichgewichtsbedingungen zwischen Abbildung 2.1: Schematische Darstellung des Phasenübergangs an einer Dampf- Flüssig-Grenzfläche im thermodynamischen Gleichgewicht. Die durch- gezogene Linie zeigt den realen, die gestrichelte Linie den idealisierten Verlauf. KAPITEL 2. THERMODYNAMISCHE GRUNDLAGEN 5 allen Phasen π über die Grenzfläche hinweg T ′ = T ′′ = ... = T π (thermisches Gleichgewicht) (2.1) p′ = p′′ = ... = pπ (mechanisches Gleichgewicht) (2.2) µ′ i = µ′′ i = ... = µπ i (stoffliches Gleichgewicht) (2.3) mit dem chemischen Potential µi für alle Komponenten i = 1, ..., n. Während ther- misches und stoffliches Gleichgewicht auch im unmittelbaren Bereich der Grenzfläche gültig sind, ist die Definition des Gleichgewichtsdrucks an der Grenzfläche nicht trivial. Für ein Dampf-Flüssig-System im thermodynamischen Gleichgewicht entspricht dieser Kernphasendruck dem Dampf- bzw. Siededruck pS des Reinstoffs oder der Mischung. Für den vereinfachten Fall einer planaren Grenzfläche in x,y-Richtung bleibt nur noch eine Richtungsabhängigkeit in der z -Koordinate. Im thermodynamischen Gleichgewicht reduziert sich der Drucktensor aufgrund des Stabilitätskriteriums (∇ ·p = 0) zu p =    pxx(z) 0 0 0 pyy(z) 0 0 0 pzz(z)    . (2.4) Es bleiben also nur noch die Hauptdiagonalelemente als Funktion der z -Koordinate übrig, alle Nebenelemente werden zu Null. Um den Richtungscharakter der Druckele- mente hervorzuheben, werden die Begriffe parallel (||) und normal (⊥) im Bezug auf die Grenzfläche eingeführt, so dass gilt p =    p||(z) 0 0 0 p||(z) 0 0 0 p⊥    . (2.5) Hierbei gilt weiterhin, dass die Normalkomponente konstant ist und damit dem Kern- phasendruck entspricht. Somit gilt die anfangs erwähnte Gleichgewichtsbedingung für das mechanische Gleichgewicht (2.2) nur für die Normalkomponente des Drucktensors. Für die Berechnung des chemischen Potential ist es zwecksmäßig, dieses auf einen Re- ferenzzustand zu beziehen. Für thermische Zustandsgleichungen als thermodynamische Modelle wird dabei der Bezug auf ein reines ideales Gas gewählt. Für einen Zustand {T, p, x} gilt µi(T, p, x) = µrein, id. Gas i (T, p0) +RT ln fi(T, p, x) p0 (2.6) mit einem willkürlich gewählten Referenzdruck p0 und der Fugazität fi. Die Fugazität 6 KAPITEL 2. THERMODYNAMISCHE GRUNDLAGEN lässt sich ansetzen über fi(T, p, x) = p ·xi ·ϕi(T, p, x) (2.7) mit dem Fugazitätskoeffizieten ϕi. Setzt man die Gleichungen (2.6) und (2.7) ineinander ein, so ergibt sich µi(T, p, x) = µrein, id. Gas i (T, p0) +RT ln p p0 ︸ ︷︷ ︸ reines ideales Gas bei T,p +RT ln xi +RT lnϕi(T, p, x) ︸ ︷︷ ︸ reale Mischung bei T,p,x . (2.8) Kombiniert man diesen Ausdruck mit Gleichung (2.3), so erhält man die Isofuagzitäts- beziehung zwischen allen Phasen π x′ iϕ ′ i = x′′ iϕ ′′ i = ... = xπ i ϕ π i . (2.9) Der Fugazitätskoeffizient lässt sich allgemein aus thermischen Zustandsgleichungen über die folgende Beziehung berechnen kT lnϕi = ∫ p p=0 ( vi − kT p ) dp . (2.10) Der Fugazitätskoeffizient entspricht daher der Abweichung des partiellen molaren Volumens vi des Stoffs i in der realen Mischung im Vergleich zum molaren Volumen der Mischung idealer Gase. Die Isofugazitätsbeziehung ist die allgemeinste Phasen- geleichgewichtsbeziehung und lässt sich universell auf die in dieser Arbeit im Fokus stehenden Phasengleichgewichte anwenden. Die Druckabhängigkeit stellt hierbei jedoch ein Problem aufgrund der Ortsabhänggkeit innerhalb der Grenzfläche dar. Deshalb ist eine alternative Formulierung im Variablenraum {T, v, x} bzw. {T, ρ} zu wählen, wie im folgenden Unterkapitel erläutert wird. Werden Grenzflächen im quasi-statischen Nicht-Gleichgewicht betrachtet, so gibt es lediglich eine über die Grenzfläche hinweg konstante Stromdichte, die Gesamtener- giestromdichte Je. Diese setzt sich zusammen aus der Wärmestromdichte und den Energiestromdichten, die über die Enthalpie an die Stoffströme gekoppelt sind. Eine detaillierte Diskussion dieser Größen wird in Kapitel 5 gegeben. 2.2 Die Helmholtzenergie und ihre Ableitungen Die in dieser Arbeit verwendete Zustandsgleichung, welche im Detail im kommenden Kapitel beschrieben wird, ist in der Helmholtzenergie F formuliert. Hierbei teilt sich die KAPITEL 2. THERMODYNAMISCHE GRUNDLAGEN 7 Helmholtzenergie in einen idealen Gasanteil und einen residuellen Anteil auf und lässt sich in Abhängigkeit der Temperatur T und des Vektors aller Komponentendichten ρ schreiben als F (T, ρ) = F ig(T, ρ) + F res(T, ρ) . (2.11) Die Berechnung des residuellen Drucks erfolgt über die thermodynamische Beziehung pres = −kT   ∂ ( F res(T,ρ) NkT ) ∂v   T,x = ρ2   ∂ ( F res(T,ρ) N ) ∂ρ   T,x (2.12) mit der Gesamtdichte ρ = ∑ i ρi. Um den Gesamtdruck zu erhalten, muss der ideale Gasanteil hinzuaddiert werden, wor- aus folgt p = pid + pres = pid + ρ2   ∂ ( F res(T,ρ) N ) ∂ρ   T,x (2.13) mit pid = ρkT . (2.14) Diese Gleichung lässt sich durch Umformung über den Kompressibilitätsfaktor Z = p ρkT umformen zu Z = 1 + ρ   ∂ ( F res(T,ρ) NkT ) ∂ρ   T,x . (2.15) Bei der Berechnung kalorischer Größen ist die Verwendung der Helmholtzenergie sehr vorteilhaft. Sie können einzig über partielle Ableitungen der Helmholtzenergie formuliert werden14. Die residuelle, molekülzahlbezogene Enthalpie ergibt sich zu ĥres(T, ρ) RT = hres(T, ρ) kT = F res(T, ρ) NkT − T   ∂ ( F res(T,ρ) NkT ) ∂T   ρ + (Z − 1). (2.16) Die isobare und isochore Wärmekapazität lassen sich ebenfalls in Abhängigkeit partieller Ableitungen der residuellen Helmholtzenergie herleiten, werden jedoch in dieser Arbeit nicht betrachtet. Das residuelle chemische Potential lässt sich aus der Fundamentalglei- chung ableiten zu µi(T, ρ) = ( ∂F ∂Ni ) T,v,Nj 6=i . (2.17) 8 KAPITEL 2. THERMODYNAMISCHE GRUNDLAGEN mitNi als Molekülzahl der Komponente i. Umgeschrieben mit der molekülzahlbezogenen Helmholtzenergie ergibt sich µres i (T, ρ) kT =   ∂ ( F res(T,ρ) V kT ) ∂ρi   T,ρj 6=i =   ∂ ( F res(T,ρ) NkT · ρ ) ∂ρi   T,ρj 6=i . (2.18) Um die Berechnungsvorschrift für den Fugazitätskoeffizienten zu erhalten, wird Glei- chung (2.8) hinzugezogen µi(T, p, x) = µid. Gas i (T, p) + kT lnϕi(T, p, x) . (2.19) Der Umstand, dass dieser Ausdruck in den Variablen (T, p, x) formuliert ist, ist nicht nachteilhaft, da für das gesamte chemische Potential gilt µi(T, p, x) = µi(T, ρ) . (2.20) Subtrahiert man folglich in Gleichung (2.19) auf beiden Seiten µid.Gas i (T, p, x) und be- rücksichtigt µid.Gas i (T, p, x) = µid.Gas i (T, v, x) + ∫ p pid vi dp = µid.Gas i (T, v, x) + ∫ p pid ( kT p ) dp , (2.21) so erhält man lnϕi(T, ρ) = µres i (T, ρ) kT − lnZ . (2.22) 9 3 PCP-SAFT Zustandsgleichung für Mischungen 3.1 Molekulares Modell Die PCP-SAFT (engl. perturbed chain polar statistical associating fluid theory) Zu- standsgleichung basiert wie einige weitere Zustandsgleichungen auf der SAFT Theorie entwickelt von Chapman et al.5,6,7. Hierbei wird ein Molekül als Kette vergröbert abstra- hiert, die aus tangential verbundenen, sphärischen Kugelsegmenten zusammengesetzt ist. Die Grundlage für die Verbindung der Kugelsegmente legte Wertheim mit seiner thermodynamischen Störungstheorie erster Ordnung (engl. thermodynamic perturbation theory of first order, TPT1)1,2,3,4. Die Kugelsegmente zwischen den Ketten interagieren untereinander durch van-der-Waals-Wechselwirkungen. Diese Wechselwirkungen werden durch ein Lennard-Jones-Potential beschrieben. Barker und Henderson28 schlagen in ih- rer Störungstheorie eine Aufteilung des Potentials in einen repulsiven Anteil φ0 ij(r) = 4εij [ (σij/r) 12 − (σij/r) 6] Θ(σij − r) (3.1) und einen attraktiven Anteil uPT ij (r) = 4εij [ (σij/r) 12 − (σij/r) 6] Θ(r − σij) (3.2) vor. Hierbei ist r der Abstand zweier Segmente i und j, σii der mittlere Segmentdurch- messer und εii der mittlere Energieparameter. Θ ist die Heaviside Stufenfunktion. Das Fluid mit dem repulsiven Potential (Gl. (3.1)) wird als Fluid aus harten Kugeln bzw. harten Kettensegmenten beschrieben, wobei der effektive Segmentdurchmesser dij(T ) einzig eine Funktion der Temperatur ist, gemäß dij(T ) = σij [ 1− 0.12 · exp ( −3εij kT )] . (3.3) Jeder Reinstoff – zunächst ohne assoziierende oder multipolare Wechselwirkungen – ist hierbei über den Durchmesserparameter σii, den Energieparameter εii sowie die Anzahl 10 KAPITEL 3. PCP-SAFT ZUSTANDSGLEICHUNG FÜR MISCHUNGEN an Segmenten mi innerhalb der Kette definiert. Für Mischungen werden die Berthelot- Lorenz Mischungsregeln für die Kreuzparameter angewendet σij = 1 2 (σii + σjj) (3.4) εij =(εiiεjj) 0,5(1− kij) (3.5) für i 6= j. Hierbei ist kij = kji ein anpassbarer Parameter, der die Segment-Segment- Wechselwirkungen ungleicher Ketten korrigiert. Die Berechnung thermodynamischer Größen erfolgt über die in Kapitel 2.2 gezeigten Ableitungen der Helmholtzenergie F . Diese lässt sich in einen idealen Gasbeitrag und einen residuellen Beitrag aufteilen F = F ig + F res . (3.6) Der residuelle Beitrag kann je nach den zu berücksichtigendenWechselwirkungsbeiträgen wiederum aufgeteilt werden in F res = FKette + F disp + F assoc + Fmultipolar . (3.7) Der Index Kette beschreiben die repulsiven Wechselwirkungen harter Ketten. Die attrak- tiven (dispersiven) Wechselwirkungen sind im Term mit dem Index disp berücksichtigt. Assoziierende (Index assoc) und multipolare (Index multipolar) Wechselwirkungen wer- den durch die letzten beiden Terme beschrieben. Auf die einzelnen Beiträge wird in den kommenden Unterkapiteln näher eingegangen. 3.2 Wechselwirkungsbeiträge 3.2.1 Hart-Kugel-Beitrag und Beitrag harter Ketten Die Wechselwirkungen von Mischungen harter Ketten werden nach der Ableitung von Chapman et al.5 formuliert. Grundlage hierfür ist die Annahme tangential irreversibel gebundener Kugeln. Die dimensionslose Helmholtzenergie für harte Ketten lautet FKette NkT = m̄ · F hs NkT − n∑ i=1 xi(mi − 1) · ln ghsii . (3.8) Hierbei ist F hs der Beitrag durch harte Kugeln, m̄ die mittlere Segmentanzahl, xi der Stoffmengenanteil einer Komponente i und ghsii die radiale Paarverteilungsfunktion für zwei Kugeln des Typs i bei Kontaktabstand in Mischungen harter Kugeln. KAPITEL 3. PCP-SAFT ZUSTANDSGLEICHUNG FÜR MISCHUNGEN 11 Boublik29 und Mansoori et al.30 entwickelten eine Theorie zur Beschreibung ei- nes Hart-Kugel-Systems resultierend im Boublik-Mansoori-Carnahan-Starling-Leland- Modell (BMCSL). Dieses beschreibt die Helmholtzenergie in dimensionsloser Form durch F hs NkT = 1 z0 · [ 3z1z2 1− z3 + z22 z3 · (1− z3)2 + ( z32 z23 − z0 ) · ln(1− z3) ] (3.9) Für harte Ketten sind die Koeffizienten zα gemittelte Werte mit zα = π 6 ρ n∑ j=1 xjmjd α j für α ∈ {0, 1, 2, 3} . (3.10) Die gemittelte Segmentanzahl m̄ errechnet sich über m̄ = n∑ i=1 ximi (3.11) als Summation über alle Komponenten n. Die radiale Paarverteilungsfunktion lässt sich für ein Hart-Kugel-System analytisch ableiten zu ghsij = 1 1− z3 + didj di + dj · 3z2 (1− z3)2 + ( didj di + dj )2 · 2z22 z3 (3.12) mit dk als temperaturabhängigem Segmentdurchmesser gemäß Gleichung (3.3). 3.2.2 Beitrag durch dispersive Wechselwirkungen Der Wechselwirkungsbeitrag durch attraktive van-der-Waals Kräfte basiert auf der Stö- rungstheorie zweiter Ordnung von Barker und Henderson. Durch die Anwendung der Ein-Fluid-Theorie kann folgende Form abgeleitet werden15 F disp NkT = −2πρI1(η̂, m̄)m2ǫσ3 − πρm̄ ( 1 + Zhc + ρ ∂Zhc ∂ρ )−1 I2(η̂, m̄)m2ǫ2σ3 (3.13) mit m2ǫσ3 = n∑ i=1 n∑ j=1 xixjmimj ( ǫij kT ) σ3 ij (3.14) und m2ǫ2σ3 = n∑ i=1 n∑ j=1 xixjmimj ( ǫij kT )2 σ3 ij . (3.15) Die Größen I1(η̂, m̄) und I2(η̂, m̄) stehen für Korrelationsintegrale, die in Abhängigkeit der mittleren Packungsdichte und der mittleren Segmentanzahl vereinfacht als Potenz- 12 KAPITEL 3. PCP-SAFT ZUSTANDSGLEICHUNG FÜR MISCHUNGEN reihen dargestellt sind. Ausformuliert lauten diese I1(η̂, m̄) = 6∑ i=0 ai(m̄)η̂i (3.16) I2(η̂, m̄) = 6∑ i=0 bi(m̄)η̂i (3.17) mit den Funktionen ai(m̄) = a0i + m̄− 1 m̄ a1i + m̄− 1 m̄ m̄− 2 m̄ a2i (3.18) bi(m̄) = b0i + m̄− 1 m̄ b1i + m̄− 1 m̄ m̄− 2 m̄ b2i . (3.19) Alle Modellkonstanten a0i, a1i, a2i sowie b0i, b1i, b2i wurden von Gross und Sadowski15 an Reinstoffdaten von n-Alkanen angepasst und als universell angesehen. Die mittlere Packungsdichte errechnet sich aus η̂ = π 6 ρ n∑ j=1 xjmjd 3 j . (3.20) Der Ausdruck ( 1 + Zhc + ρ∂Zhc ∂ρ ) lässt sich umschreiben zu ( 1 + Zhc + ρ ∂Zhc ∂ρ ) = ( 1 + m̄ 8η̂ − 2η̂2 (1− η̂)4 + (1− m̄) 20η̂ − 27η̂2 + 12η̂3 − 2η̂4 [(1− η̂)(2− η̂)]2 ) . (3.21) Um einen weiteren Freiheitsgrad in Form eines zusätzlichen anpassbaren Parameters lij zu bekommen, haben Tang und Gross31 den Term m2ǫσ3 in Gleichung (3.13) erweitert m2ǫσ3 = n∑ i=1 n∑ j=1 xixjmimj ( ǫij kT ) σ3 ij+ n∑ i=1 ximi ( n∑ j=1 xjmjσij [( ǫij kT ) lij ]1/3 )3 . (3.22) 3.2.3 Beitrag durch assoziierende Wechselwirkungen Chapman et al.32 entwickelten auf Basis der Arbeiten von Wertheim1,2,3,4 einen Helm- holtzenergiebeitrag durch assoziierende Wechselwirkungen F assoc RT = n∑ i=1 xi [ ∑ Ai ( lnXAi − XAi 2 ) + 1 2 Mi ] . (3.23) Die Modellbetrachtung hinter diesem Ansatz ist, dass ein Monomer-Kugelsegment oder eine Kette mit mehreren Kugelsegmenten mehrere Bindungsstellen A, B, C etc. aufwei- sen kann. In der in dieser Arbeit angewendeten Form wurde die Anzahl der Bindungs- KAPITEL 3. PCP-SAFT ZUSTANDSGLEICHUNG FÜR MISCHUNGEN 13 stellen auf zwei, also A und B beschränkt. In Gleichung (3.23) beschreibt die Größe XAi den Stoffmengenanteil der Komponente i, der nicht an der Bindungsstelle A gebunden ist. Berechnen lässt sich dies über XAi =  1 +NAv n∑ j=1 ∑ Bj ρjX Bj∆AiBi   −1 . (3.24) Hierbei beschreibt ∑ Bj die Summe über alle Bindungsstellen an Molekül j, XBj den Stoffmengenanteil der Komponente j, der nicht an der Bindungsstelle B gebunden ist, und ∆AiBi die Assoziierungsstärke. Diese ist definiert über ∆AiBi = d3ijg hs ij κ AiBj [ exp ( εAiBj/kT ) − 1 ] (3.25) mit zwei weiteren Reinstoffparametern κAiBj als Assoziierungsvolumen und εAiBj als Assoziierungsenergie. Die Definition des mittleren Segmentdurchmessers dij entspricht Gleichung (3.3). Die radiale Paarverteilungsfunktion harter Kugeln ghsij ist durch Glei- chung (3.12) beschrieben. 3.2.4 Beitrag durch multipolare Wechselwirkungen Um multipolaren Wechselwirkungen zwischen Molekülketten Rechnung zu tragen, wer- den innerhalb der Zustandgleichung Beiträge durch dipolare und quadrupolare Anteile berücksichtigt. Beide Anteile werden in Form einer Störungstheorie dritter Ordnung beschrieben. Für den dipolaren Anteil der Helmholtzenergie entwickelten Gross und Vrabec17 ein Modell gemäß F dd,EoS NkT = F dd,EoS 2 /(NkT ) 1− F dd,EoS 3 /F dd,EoS 2 (3.26) mit dem Störungsterm zweiter Ordnung F dd,EoS 2 NkT = −πρ n∑ i=1 n∑ j=1 xixj εii kT εjj kT σiiσjj σij nµ,inµ,jµ ∗ i 2µ∗ j 2JDD 2,ij (3.27) und dem Störungsterm dritter Ordnung F dd,EoS 3 NkT = −4π2 3 ρ2 n∑ i=1 n∑ j=1 n∑ k=1 xixjxk εii kT εjj kT εkk kT σ3 iiσ 3 jjσ 3 kk σijσikσjk nµ,inµ,jnµ,kµ ∗ i 2µ∗ j 2µ∗ k 2JDD 3,ijk . (3.28) Hierbei steht nµ für die Anzahl an Dipolmomenten pro Molekül, welche jedoch in dieser Arbeit immer gleich 1 gesetzt ist. Das Quadrat des dimensionslosen Dipolmoments einer 14 KAPITEL 3. PCP-SAFT ZUSTANDSGLEICHUNG FÜR MISCHUNGEN beliebigen Komponente l ist definiert mit µ∗ l 2 = µ2 l mlεllσ3 ll . (3.29) JDD 2,ij und JDD 3,ijk beschreiben Integrale über die Paarverteilungsfunktion bzw. die Drei- Körper-Korrelationsfunktion des Referenzfluids und werden analog zur Vorgehensweise im Beitrag durch dispersive Wechselwirkungen durch Potenzreihen dargestellt mit JDD 2,ij = 4∑ n=0 ( an,ij + bn,ij εij kT ) ηn (3.30) und JDD 3,ijk = 4∑ n=0 cn,ijkη n . (3.31) Die Koeffizienten berechnen sich aus an,ij =a0n + mij − 1 mij a1n + mij − 1 mij · mij − 2 mij a2n (3.32) bn,ij =b0n + mij − 1 mij b1n + mij − 1 mij · mij − 2 mij b2n (3.33) cn,ijk =c0n + mijk − 1 mijk c1n + mijk − 1 mijk · mijk − 2 mijk c2n (3.34) mit den gemittelten Segmentanzahlen mij = min ( 2, (mimj) 1/2 ) (3.35) und mijk = min ( 2, (mimjmk) 1/3 ) . (3.36) Die Einschränkung der gemittelten Segmentanzahlen auf 2 stellt dabei sicher, dass sich das Dipolmoment nur auf maximal zwei Segmente ausdehnt. Grundlage für den Beitrag der Helmholtzenergie durch dipolare Wechselwirkungen war die Formulierung des Beitrags der Helmholtzenergie durch quadrupolare Wechselwir- kungen16 gemäß] FQQ NkT = FQQ 2 /(NkT ) 1− FQQ 3 /FQQ 2 (3.37) mit dem Störungsterm zweiter Ordnung FQQ 2 NkT = −π ( 3 4 )2 ρ n∑ i=1 n∑ j=1 xixj εii kT εjj kT σ5 iiσ 5 jj σ7 ij nQ,inQ,jQ ∗ i 2Q∗ j 2JQQ 2,ij (3.38) KAPITEL 3. PCP-SAFT ZUSTANDSGLEICHUNG FÜR MISCHUNGEN 15 und dem Störungsterm dritter Ordnung FQQ 3 NkT = π 3 ( 3 4 )2 ρ n∑ i=1 n∑ j=1 xixj ( εii kT )3/2 ( εjj kT )3/2 σ 15/2 ii σ 15/2 jj σ12 ij × nQ,inQ,jQ ∗ i 3Q∗ j 3JQQ 3,ij + 4π2 3 ( 3 4 )2 ρ2 n∑ i=1 n∑ j=1 n∑ k=1 xixjxk × εii kT εjj kT εkk kT σ5 iiσ 5 jjσ 5 kk σ3 ijσ 3 ikσ 3 jk nQQ,inQQ,jnQ,kQQ∗ i 2QQ∗ j 2QQ∗ k 2JQQ 3,ijk . (3.39) Die Integrale JQQ 2,ij über die Paarverteilungsfunktion und JQQ 3,ijk über die Drei-Körper- Korrelationfunktion werden analog zu den Gleichungen (3.30) und (3.31) formuliert JQQ 2,ij = 4∑ n=0 ( an,ij + bn,ij εij kT ) ηn (3.40) JQQ 3,ijk = 4∑ n=0 cn,ijkη n . (3.41) Der Ansatz für die Koeffizienten erfolgt ebenfalls analog zu den Gleichungen (3.32) bis (3.34). Um die Anzahl an Koeffizienten nicht weiter zu erhöhen, wird das Integral JQQ 3,ij = 0 gesetzt. Die Anzahl der Quadrupolmomente pro Molekül nQQ wird in dieser Arbeit gleich 1 gesetzt. Analog zu den gezeigten Beiträgen zu quadupolaren und dipolaren Wechselwirkungen ist für Mischungen aus dipolaren und quadrupolaren Komponenten ebenfalls der Beitrag zu dipolar-quadrupolaren Wechselwirkungen (DQ) zu beschreiben. Der Übersichtlichkeit halber wird auf die vollständige Formulierung der Gleichungen an dieser Stelle verzichtet. 16 4 Dichtefunktionaltheorie 4.1 Das Großkanonische Potential Ω Zur Untersuchung von Grenzflächeneigenschaften ist ein thermodynamisches Potential mit den natürlichen Variablen chemisches Potential µ und Temperatur T vorzuziehen. Da diese beiden Variablen im Gleichgewichtszustand eines Systems über die Grenzfläche hinweg konstant sind, lassen sich daraus Berechnungen stark vereinfachen. Das Poten- tial mit diesen natürlichen Variablen ist das Großkanonische Potential Ω(T, V, µ). Aus der Fundamentalgleichung für die innere Energie U(S, V,N) lässt dich über Legendre- Transformation das Großkanonische Potential in Abwesenheit eines externen Feldes de- finieren als Ω(T, V, µ) = F (T, V,N)− µN (4.1) mit der Helmholtzenergie F (T, V,N) und der Gesamtmolekülzahl N. Volumen und Gesamtmolekülzahl lassen sich für ein inhomogenes System beschreiben als V = ∫ dr (4.2) N = ∫ ρ(r) dr (4.3) mit r als Raumvektor über alle drei Raumkoordinaten. In Funktionalschreibweise ergibt sich damit für das Großkanonische Potential einer Multikomponentenmischung Ω[ρk(r̃);T, µ] = F [ρk(r̃)]− n∑ i=1 µi ∫ dr̃ ρi(r̃) . (4.4) Hierbei ist der Index k ein generischer Index, der für alle Komponenten n steht (k = 1, . . . , n). Die Größe ρi(r̃) steht für die Molekulardichte mit r̃ als Vektor, der sowohl die Position als auch die Konfiguration von Molekülen beschreibt. Befindet sich ein System im thermodynamischen Gleichgewicht, so wird das Großkano- nische Potential minimal 0 = δF [ρk] δρi(r̃) − µi ∀i . (4.5) KAPITEL 4. DICHTEFUNKTIONALTHEORIE 17 Diese Gleichungen stellen im Allgemeinen die Berechnungsvorschrift für den Verlauf der Dichten aller Stoffe über die Grenzfläche hinweg dar. Im Abschnitt 4.2.1 wird eine Modi- fikation dieser Gleichung gezeigt, die die Berechnung des Dichteprofils weiter vereinfacht und zugänglich macht. Die Berechnung der Grenzflächenspannung erfolgt direkt aus der Differenz des Großka- nonischen Potentials Ω [ρk] in Abhängigkeit der Gleichgewichtsdichteprofile ρk(r̃) und des Großkanonischen Potentials in der Kernphase. Die Berechnunggleichung lautet γ = 1 A (Ω [ρk(r̃)]− Ωbulk) (4.6) mit der Fläche der planaren Grenzfläche A. Die Gleichung (4.6) kann mit Gleichung (4.4) umgeschrieben werden. Für eine planare Grenzfläche ergibt dies γDFT = ∫ { f [ρk(r̃)]− n∑ i=1 µiρi(r̃) + pbulk } dr̃ (4.7) mit der Helmholtzenergiedichte f [ρk(r̃)] definiert über F [ρk(r̃)] = ∫ f [ρk(r̃)] dr̃. Berücksichtigt man zusätzlich den Effekt durch Kapillarwellen wird ein zusätzlicher Faktor multiplikativ hinzugefügt27,33,34. Dies ergibt die makroskopische Grenzflächen- spannung γ = γDFT ( 1 + 3 8π T Tc 1 (2.55)2 1 κ )−1 (4.8) mit dem als konstant angenommenen Amplitudenverhältnis der Kapillarwellen κ = 0, 24 4.2 Funktionalbeiträge zur Helmholtzenergie Die Eigenschaft, die Helmholtzenergie F aus verschiedenen Beiträgen additiv zusammen zu setzen, wurde in Abschnitt 3.2 diskutiert. Diese Additivität wird auch in der Funk- tionalbetrachtung innerhalb des DFT Formalismus beibehalten. Folgende Aufteilung in idealen Gasbeitrag und residuellen Anteil kommt ganz allgemein zur Anwendung F [ρk(r̃)] = F ig[ρk(r̃)] + F res[ρk(r̃)] . (4.9) Der residuelle Beitrag lässt sich weiterhin in die Beiträge durch Wechselwirkungen harter Kugeln (eng. hard spheres, Index hs), Kettenformation (Index Kette), dispersive Wech- selwirkungen (Index disp), assoziierende Wechselwirkungen (Index assoc) sowie polare 18 KAPITEL 4. DICHTEFUNKTIONALTHEORIE Wechselwirkungen (Index polar) aufteilen F [ρk(r̃)] = F ig[ρk(r̃)]+F hs[ρk(r̃)]+FKette[ρk(r̃)]+F disp[ρk(r̃)]+F assoc[ρk(r̃)]+F polar[ρk(r̃)] (4.10) Analog ergibt sich das chemische Potential aus den additiven Beiträgen µ = µig + µhs + µKette + µdisp + µassoc + µpolar . (4.11) Auf die einzelnen Beiträge und die getroffenen Vereinfachungen hierzu wird in den fol- genden Abschnitten näher eingegangen. 4.2.1 Idealer Gasbeitrag Die Funktionalbeschreibung des idealen Gasbeitrags der Helmholtzenergie in dimensi- onsloser Form ist folgendermaßen definiert F ig[ρk(r)]/kT = ∫ n∑ i=1 ρi(r) { ln [ ρi(r)Λ 3 i ] − 1 } dr (4.12) mit Λ als DeBroglie-Wellenlänge. Hierbei wurde der Übergang von ρk(r̃) auf ρk(r) über die Gleichung ρk(r̃) ≈ ρk(r)fk(ω) (4.13) vorgenommen, also das Produkt aus molekularer Dichte in Abhängigkeit des Ortes und einer Orientierungs- und Konformationsverteilungsfunktion fk(ω), die als orts- und dichteunabhängig angenommen wird. Alle Freiheitsgrade durch Orientierung und Kon- formation werden in der DeBroglie-Wellenlänge berücksichtigt35, so dass der Übergang wie gezeigt vollzogen werden kann. Die Funktionalableitung des idealen Gasbeitrags in Abhängigkeit des Dichteprofils für eine Komponente i ergibt δF ig[ρk(r)]/kT δρi(r) = ln [ ρi(r)Λ 3 i ] . (4.14) Die Kombination der Gleichungen (4.14), (4.5) und (4.11) liefert die folgende Iterations- vorschrift für die Dichteiteration einer Komponente i, wie sie für alle DFT Berechnungen in dieser Arbeit genutzt wurde, ρi(r) = ρli · exp ( µres,l i /kT − δF res[ρk]/kT δρi(r) ) ∀i . (4.15) KAPITEL 4. DICHTEFUNKTIONALTHEORIE 19 Als Referenz des chemischen Potentials wurde hierzu die Flüssigphase l gewählt. Die Wahl ist durch die Ortsunabhängigkeit des chemischen Potentials über die Grenzfläche im thermodynamischen Gleichgewicht keiner Einschränkung unterlegen. 4.2.2 Hart-Kugel-Beitrag - Fundamental Measure Theorie Eine genaue Beschreibung der Helmholtzenergie von Mischungen aus harten Kugeln wurde durch die Fundamental Measure Theorie (FMT) von Rosenfeld36 möglich. Die durch Roth et al.37 und Yu und Wu38 unabhängig voneinander entwickelte modifizierte FMT macht die Theorie für inhomogene Mischungen harter Kugeln mit der Boublik- Mansoori-Carnahan-Starling-Leland (BMCSL) Theorie29,30 kompatibel. Der Ausdruck für den Helmholtzenergiebeitrag einer Mischung harter Kugeln in dimensionsloser Form ist F hs[ρ] kT = ∫ fhs (n) kT dr . (4.16) Hierbei ist die Helmholtzenergiedichte fhs(n) nur eine Funktion von vier Dichten nα mit α ∈ {0, . . . , 3}, die wiederum Funktionale der Dichte sind. Weiterhin gehen zwei gewichtete Vektordichten nα mit α ∈ {V 1, V 2} ein. Der gesamte Ausdruck ergibt sich zu fhs(n) kT = − n0 ln(1− n3) + n1n2 − nV 1 ·nV 2 1− n3 + ( n3 2 − 3n2nV 2 ·nV 2 ) n3 + (1− n3) 2 ln(1− n3) 36πn2 3(1− n3)2 . (4.17) Die Definition der Dichten nα und nα sind von Gross27 für planare Grenzflächen gege- ben. Die Funktionalableitung des Energiebeitrags für eine Komponente i ergibt für eine pla- nare Grenzfläche ausgewertet δF hs[ρk]/kT δρi(x) = mi ∫ di/2 −di/2 [ 1 di ∂φ ∂n0 + 1 2 ∂φ ∂n1 + πdi ∂φ ∂n2 + π ∂φ ∂n3 ( (di/2) 2 − z′2 ) ] dz′ +miez ∫ di/2 −di/2 [ 1 di ∂φ ∂nV 1 z′ + 2π ∂φ ∂nV 2 z′ ] dz′ . (4.18) Hierbei ist φ = fhs/kT . Für den Durchmesser di wird der temperaturabhängige, effek- tive Hart-Kugel-Durchmesser nach Gleichung (3.3) verwendet. Die in Gleichung (4.18) vorkommenden Ableitungen ∂φ/∂nα und ∂φ/∂nβ können für planare Grenzflächen analog zu den Ausführungen von Gross27 für jede Komponente ausgewertet werden. Dies gewährleistet die Konsistenz zur PCP-SAFT Zustandsgleichung. 20 KAPITEL 4. DICHTEFUNKTIONALTHEORIE 4.2.3 Kettenbeitrag Tripathi und Chapman39,40 entwickelten die funktionale Form des Helmholtzenergie- beitrags durch Kettenbildung nicht-sphärischer Kugeln. Diese basiert auf Wertheims Störungstheorie erster Ordnung1,2,3,4. Der große Vorteil der Formulierung von Tripathi und Chapman ist die Wahl des Referenzzustandes. Dieser beschreibt die Mischung nicht- gebundener sphärischer Segmente im idealen Gaszustand. Prinzipiell müssten hierzu alle Bedingungen der Bindungslängen und Winkelverteilungen der beteiligten Teilchen im idealen Gaszustand beschrieben werden. Indem die Wertheim’sche Theorie von einem Fluid ungebundener sphärischer Kugeln – hier im Kontext von Hartkugeln – ausgeht, ist es nicht erforderlich, Details zu intramolekularen Bindungen zu definieren. Eine de- taillierte Aufarbeitung und Erklärung der Wertheim’schen Theorie wurde von Zmpitas und Gross41 publiziert. Im hier vorgestellten DFT Formalismus wird der Fall von gemittelten, einheitlichen Segmenten in einer Kettenstruktur betrachtet. Diese Mittelung verhindert, dass eine Unterscheidung gemacht werden muss zwischen Dichteprofilen von Segmenten, die sich an den Enden einer Kette befinden, im Vergleich zu Segmenten innerhalb der Kette. Alle Dichteprofile ρi(r) aller Segmenttypen einer homogenen Kette i werden als gleichwertig behandelt. Damit ergibt sich die einfache Form des Kettenterms FKette[ρk(r)]/kT = n∑ i=1 (mi − 1) ∫ ρi(r) {ln (ρi(r))− 1} dr − n∑ i=1 (mi − 1) ∫ ρi(r) { ln [ yddii (ρ̄k(r)) λi(r) ] − 1 } dr . (4.19) Die Form der Gleichung mit zwei Termen spiegelt den Anteil im idealen Gaszustand und den residuellen Anteil sichtbar wieder. yddii (ρ̄k(r)) ist die Kavitätskorrelationsfunktion, welche am Berührungspunkt zweier nicht gebundener Segmente di ausgewertet wird. Im Berührungspunkt entspricht der Wert von yddii dem der radialen Paarverteilungsfunktion. Ein starkes Modell für yddii (r− r′) zu entwickeln ist nicht trivial, so dass der Ansatz von Kierlik und Rosinberg42 sowie von Tripathi und Chapman39 angewendet wird. Hierbei wird yddii (r− r′) angenähert durch yddii (r− r′) ≈ yddii (ρ̄k(r)) (4.20) mit der mittleren Dichte ρ̄k(r). Wertet man die Boublik-Mansoori-Carnahan-Starling- Leland (BMCSL) Theorie29,30 für eine ortsfeste, mittlere Dichte ρ̄k(r) aus, so ergibt sich yddii (ρ̄k(r)) = 1 1− ζ̄3 + 3/2 diζ̄2 (1− ζ̄3)2 + 0.5(diζ̄2) 2 (1− ζ̄3)3 (4.21) KAPITEL 4. DICHTEFUNKTIONALTHEORIE 21 mit gemittelten Packungsdichten (n = (2, 3)) ζ̄n(r) = π 6 n∑ i=1 ρ̄i(r)mid n i . (4.22) ρ̄i(r) ist die mittlere Dichte im ausgeschlossenen Volumen eines Segments and einer festen Position r als ρ̄i(r) = 3 4πd3i ∫ ρi(r ′) Θ(di − |r− r′|) dr′ . (4.23) Die Größe λi(r) in Gleichung (4.19) beschreibt die gemittelte Dichte am Berührungs- punkt um ein Segment der Kette i am Ort r über λi(r) = 1 4πd2i ∫ ρi(r ′) δ(di − |r− r′|) dr (4.24) mit der Dirac-Funktion δ. Die Funktionalableitung des Kettenterms ergibt sich zu δFKette[ρk]/kT δρi(r) =(mi − 1) ln (ρi(r)))− (mi − 1) { ln [ yddii (ρ̄k(r))λi(r) ] − 1 } − n∑ j=1 (mj − 1) ∫ ρj(r ′) ∂ ln yddjj (ρ̄k(r ′)) ∂ρ̄i(r′) 3 4πd3i Θ(di − |r′ − r|) dr − (mi − 1) ∫ ρi(r ′) 1 λi(r′) 1 4πd2i δ(di − |r′ − r|) dr . (4.25) 4.2.4 Entwicklung eines neuen Helmholtzenergiefunktionals dispersiver Wechselwirkungen für Mischungen Bisherige Studien27 beschäftigten sich nur mit der Entwicklung eines Helmholtzener- giefunktionals dispersiver Wechselwirkungen für Reinstoffe. Für die Entwicklung eines Helmholtzenergiebeitrags für Mischungen dient der in allgemeiner Form formulierte, stö- rungstheoretische Ansatz F 1PT[ρk]/kT = 1 2 n∑ i=1 n∑ j=1 mi∑ αi=1 mj∑ βj=1 ∫ ∫ ρi,αi (r) ρj,βj (r′) × ghciαijβj (r, r′) φiαijβj (r, r′) kT drdr′ + . . . (4.26) Betrachtet man nur den (dargestellten) ersten Term und vernachlässigt alle weiteren Terme, so ergibt sich eine Störungstheorie erster Ordnung. Der gesamte Helmholtzener- 22 KAPITEL 4. DICHTEFUNKTIONALTHEORIE giebeitrag ergibt sich aus den Summationen über alle Komponenten i und j sowie über die Summationen über alle Segmente αi einer Kette mit mi Segmenten bzw. βj einer Kette mit mj Segmenten. ghciαijβj (r, r′) beschreibt die radiale Paarverteilungsfunktion zwischen der Segmentkombination αi ↔ βj eines Referenzfluids harter Ketten an jeder Paaranordnung im Raum, welche das Wechselwikungspotential φiαijβj (r, r′) zwischen den Segmenten gewichtet. r und r′ beschreiben die Positionsvektoren der jeweiligen Seg- mente. In dieser allgemeinen Form ist eine analytische Auswertung nicht möglich. Daher werden drei Vereinfachungen getroffen. Die einzelnen Segmente innerhalb einer Kette sollen als ununterscheidbar im Bezug auf ihre Ortsverteilung gelten, so dass sich die Segmentdich- ten zu ρi,α(r) = ρi(r) (4.27) ρj,βj (r′) = ρj(r ′) (4.28) reduzieren. Des Weiteren werden Moleküle als Segmentketten aufgefasst, die aus gleich- artigen Segmenten bestehen, so dass für die Wechselwirkungspotentiale gilt φiαijβj (r, r′) = φi,j(r, r ′) . (4.29) Mit den oben genannten Annahmen bezieht sich die Doppelsumme ∑mi αi=1 ∑mj βj=1 über Segment-Paarungen α & β in Gleichung (4.26) nur noch auf die radiale Paarverteilungs- funktion, so dass sich die Definition einer gemittelten Größe anbietet, gemäß 1 mimj mi∑ αi=1 mj∑ βj=1 ghciαijβj (r, r′) ≡ ghci,j (r, r ′) . (4.30) Um die numerische Auswertung der radialen Paarverteilungsfunktion weiter zu verein- fachen, wird diese über die Form einer Drei-Fluid-Theorie angenähert mit ghci,j (r, r ′) ≈ ghc (r̂, η̂, mij) , (4.31) wobei ein inhomogenes Fluid über die Mittelung der Packungsdichte an die eines ho- mogenen Fluids approximiert wird43. Die Abhängigkeit besteht somit aus dem Abstand r̂ = |r− r′| zwischen den Positionsvektoren, der gemittelten Packungsdichte η̂ sowie der Segmentanzahl mij. Fasst man alle beschriebenen Annahmen zusammen, so folgt für den Helmholtzenergiebeitrag F 1PT[ρk]/kT = 1 2 n∑ i=1 n∑ j=1 mimj ∫ ∫ ρi(r) ρj(r ′) ghc (r̂, η̂, mij) φPT ij (r̂) kT drdr′ . (4.32) KAPITEL 4. DICHTEFUNKTIONALTHEORIE 23 Für eine beliebige Kombination aus zwei Komponenten i-j ergibt sich das Korrelations- integral damit aus drei Fluiden. Jedes Fluid besitzt seine eigene Segmentanzahl mii, mjj und mij, wobei jedoch alle drei Fluide dieselbe gemittelte Packungsdichte η̂ haben. Die Mittelung erfolgt über die Packungsdichte bei r und r′ durch η̂ = 1 2 (η(r) + η(r′)) = 1 2 π 6 n∑ i=1 mid 3 i (ρi(r) + ρi(r ′)) . (4.33) Die Segmentanzahl wird arithmetisch zwischen den Segmentanzahlen zweier reiner Flui- de i und j gemittelt mij = 1 2 (mii +mjj) . (4.34) Die Paarverteilungsfunktion eines Hart-Kugel Fluids wird durch die Ornstein-Zernicke Gleichung mit der Percus-Yevick Näherung abgeschätzt. Der analytische Ausdruck hier- zu wurde von Tang und Lu44 basierend auf der Arbeit von Chiew45 entwickelt. Der Ansatz einer Drei-Fluid-Theorie eignet sich im Kontext der Dichtefunktionalberech- nung besonders, da die radiale Paarverteilungsfunktion für die drei Kettenlängen vorab für verschiedene Abstände r̂ und Packungsdichten η̂ berechnet werden kann. Durch ein- fache Interpolation zwischen den vorab ermittelten Werten kann der Rechenaufwand in jeder Iterationsschleife stark verringert werden. Bei einer Ein-Fluid-Theorie müsste die Paarverteilungsfunktion ghc im Gegensatz dazu für jede Zusammensetzung entlang der Grenzfläche einzeln berechnet werden. Reduziert man Gleichung (4.32) auf ein homogenes Fluid, so fällt das entwickelte Funk- tional nicht auf die Helmholtzenergie des dispersiven Wechselwirkungsanteils aus der PCP-SAFT Zustandsgleichung zurück. Die Differenz zwischen dem Wechselwirkungsan- satz aus der PCP-SAFT Zustandsgleichung und der Funktionalformulierung ist jedoch recht gering, so dass sie durch eine lokale Dichteapproximation, wie von Gloor et al.46,47 vorgeschlagen, angewandt werden kann. Der volle Dispersionsterm formuliert in der Helmholtzenergiedichte ergibt sich zu fdisp[ρk(r)] = f 1PT[ρk(r)] + (fdisp,PC-SAFT(ρk)− f 1PT(ρk)) . (4.35) Die eingeklammerte Differenz auf der rechten Seite der Gleichung (4.35) dient dazu, die Konsistenz mit der PCP-SAFT Zustandsgleichung für Kernphasen zu erreichen. Im eingeklammerten Term kommt die lokale Dichteapproximation (LDA) zum Einsatz, die für den gesamten Dispersionsterm eigentlich eine zu grobe Näherung darstellt. Für die relativ kleine Korrektur ist die LDA jedoch ausreichend. Nach der LDA wird die Helmholtzenergiedichte einfach für die lokalen Werte des Dichteprofils ρk = ρk(r) so ausgewertet, als läge dort ein homogenes Fluid mit der entsprechenden Dichte vor. 24 KAPITEL 4. DICHTEFUNKTIONALTHEORIE Die Funktionalableitung bezogen auf die Dichte einer Komponente i ist δF disp[ρk]/kT δρi(r) = δF 1PT[ρk]/kT δρi(r) + ( µdisp,PC-SAFT i kT − µ1PT i kT ) (4.36) mit δF 1PT[ρk]/kT δρi(r) = n∑ j=1 mimj ∫ ρj(r ′) ghc (r̂, η̂, mij) φPT ij (r̂) kT dr′ + 1 2 n∑ j=1 n∑ k=1 mjmk ρk(r) ∫ ρj(r ′) ∂ghc(r̂, η̂, mjk) ∂η̂ π 6 mid 3 i φPT jk (r̂) kT dr′ (4.37) Die Integration in Gleichung (4.37) wird numerisch bis zu einem Abschneideradius rc = 9σ ausgeführt. Außerhalb dieses Abschneideradius kommen Tail-Korrekturen zur Anwendung27. 4.2.5 Entwicklung eines neuen Helmholtzenergiefunktionals dipolarer Wechselwirkungen In vorangegangenen Studien27 wurde für stark dipolare Stoffe systematische Abweichun- gen konstatiert. Es hat sich die These gefestigt, dies sei auf die Vernachlässigung der Orientierungsverteilung zurück zu führen. In den bestehenden Funktionalen wird ein iso- tropes Verhalten unterstellt. Ziel war es daher, anisotrope Effekte im DFT Formalismus abbilden zu können. Dazu muss Gleichung (4.4) um winkelabhängige Beiträge erweitert werden. Für ein Reinstoffsystem lässt sich das Großkanonische Potential formulieren als Ω = F [ρ̂(r,ω)]− ∫ ∫ µ̂(r,ω)ρ̂(r,ω) dωdr . (4.38) Die Größe ρ̂(r,ω) beschreibt sowohl die Position eines Teilchens mit dem Ortsvektor r als auch die Ausrichtung des Dipolmoments ω. Im Vergleich zur isotropen Formulierung des Großkanonische Potentials ist hier auch das chemische Potential winkelabhängig. Dies zeigt die Analyse des Systems im thermodynamischen Gleichgewicht, bei dem wie im isotropen Fall das Großkanonische Potential minimal wird. Es gilt δΩ δρ̂(r,ω) ∣ ∣ ∣ ∣ eq = 0 . (4.39) KAPITEL 4. DICHTEFUNKTIONALTHEORIE 25 Die Funktionalableitung definiert also auch in der erweiterten Form das chemische Po- tential δF [ρ̂(r,ω)] δρ̂(r,ω) ∣ ∣ ∣ ∣ eq = µ̂(r,ω) . (4.40) Die Dichte ρ̂(r,ω) kann als Produkt zerlegt werden ρ̂(r,ω) = ρ(r)α(r,ω) , (4.41) mit der bekannten Molekülzahldichte ρ(r) = ∫ ρ̂(r,ω) dω , (4.42) welche durch Integration von ρ̂(r,ω) über alle Winkel erhalten wird. α(r,ω) bezeich- net die Orientierungsverteilungsfunktion des Dipolmoments. Die Zerlegung in Gleichung (4.41) ist naheliegend, weil einerseits die gewöhnliche Molekülzahldichte ρ(r) genutzt wird und andererseits die Orientierungsverteilungsfunktion eine sehr geschickte Normie- rung hat. Dies zeigt sich aus der Kombination aus den Gleichungen (4.41) und (4.42), woraus folgt ∫ α(r,ω) dω = 1 . (4.43) Für anisotrope Fluide reduziert sich die Orientierungsverteilungsfunktion zu einer Kon- stanten, αiso(r,ω) = 1 4π . Die erweiterte Dichtefunktion ergibt sich im isotropen Fall demnach zu ρ̂iso(r,ω) = 1 4π ρ(r) . (4.44) Die Analyse nicht winkelabhängiger Helmholtzenergiebeiträge ist im Anhang A in die- sem Kontext diskutiert. In Anlehnung an die Formulierung des Helmholtzenergiefunktionals zur Beschreibung dispersiver Wechselwirkungen erfolgt auch der Ansatz zur Beschreibung dipolarer Wech- selwirkungen auf Basis einer Störungstheorie erster Ordnung F dd[ρ̂(r,ω)] = 1 2 ∫ 1 0 ∫∫∫∫ ρ̂(r,ω)ρ̂(r′,ω′)gλ(r, r ′,ω,ω′) × uDD(r, r ′,ω,ω′) drdr′dωdω′dλ (4.45) Für sphärische Moleküle mit statischen Punktdipolen und harter Repulsion ist dieses störungstheoretische Funktional exakt. Zur praktischen Anwendung sind Näherungen notwendig. Die dipolaren Wechselwirkungen werden durch das Potential uDD(r, r ′,ω,ω′) beschrie- ben. Im folgenden Unterkapitel wird näher auf das verwendete Wechselwirkungspotential eingegangen. Weiterhin bezeichnet gλ die Paarverteilungsfunktion für ein Wechselwir- 26 KAPITEL 4. DICHTEFUNKTIONALTHEORIE kungspotential, das durch den Wert λ gekennzeichnet ist. Diese Paarverteilungsfunktion ist unbekannt und muss näherungsweise modelliert werden. Hier wird ein Ansatz von Frodl und Dietrich48 erweitert gλ(r, r ′,ω,ω′) ≈ ghs(r, r′) · exp(−βλuDD(r, r ′,ω,ω′)) . (4.46) Hierbei ist ghs(r, r′) die Paarverteilungsfunktion eines Hartkugelsystems als Referenz. Die Integration über den Kopplungsparameter λ liefert F dd[ρ̂(r,ω)] = 1 2β ∫∫∫∫ ρ̂(r,ω)ρ̂(r′,ω′) × ghs(r, r′) · ( 1− exp(−βuDD(r, r ′,ω,ω′)) ) drdr′dωdω′ . (4.47) Dieser Gleichung liegt zunächst die Näherung zu Grunde, dass die Kavitätskorrelati- onsfunktion für eine gemittelte Dichte eines homogenes Fluids angenähert werden kann mit y(r, r′,ω,ω′) = y(r̂, ρ̄,ω,ω′) . (4.48) Weiterhin wird eine Entwicklung erster Ordnung in der Kavitätskorrelationsfunktion y(r̂, ρ̄,ω,ω′) um das Hartkugelfluid gemacht mit y(r̂, ρ̄,ω,ω′) = yhs(r̂, ρ̄) . (4.49) Die Reihenentwicklung in y hat bessere Konvergenzeigenschaften als die Entwicklung in der radialen Paarverteilungsfunktion g(r̂, ρ̄,ω,ω′), was auch daran zu erkennen ist, dass der erste Ordnungsterm der Entwicklung in g(r̂, ρ̄,ω,ω′) um das Hartkugelfluid Null ergibt. Durch die hohe Anzahl an durchzuführenden Winkelintegrationen in Gleichung (4.47) ist es mit vertretbarem Rechenaufwand nicht möglich, das Integral ohne analytische Aus- wertung einiger Integrale auszuführen. Dies ist jedoch für Integranden unmöglich, die das Wechselwirkungspotential uDD beinhalten. Daher sind Vereinfachungen und Nähe- rungen erforderlich. Der Exponentialausdruck in Gleichung (4.47) wird durch die Reihe exp(−x) = 1− x+ 1 2 x2 − 1 6 x3 + ... (4.50) entwickelt und ergibt nach Integration über λ F dd[ρ̂(r,ω)] = 1 2 ∫∫∫∫ ρ̂(r,ω)ρ̂(r′,ω′)ghs(r, r′)UDD(r, r ′,ω,ω′) drdr′dωdω′ . (4.51) KAPITEL 4. DICHTEFUNKTIONALTHEORIE 27 Zur Verbesserung der Übersichtlichkeit werden die Terme des Wechselwirkungspotentials zusammengefasst UDD(r, r ′,ω,ω′) = uDD(r, r ′,ω,ω′)− 1 2 βu2 DD(r, r ′,ω,ω′) + 1 6 β2u3 DD(r, r ′,ω,ω′) . (4.52) Die Funktionalableitung von Gleichung (4.51) nach der Dichte ρ̂(r,ω) liefert δF dd[ρ̂(r,ω)] δρ̂(r,ω) = ∫ ∫ ρ̂(r′,ω′) ( ghs(r, r′) + 1 2 ρ(r) ∂ghs(r, r′) ∂ρ̄(r, r′) ) UDD(r, r ′,ω,ω′) dr′dω′ (4.53) mit der mittleren Dichte ρ̄ = 1 2 (ρ(r) + ρ(r′)). Wie im Falle des Helmholtzenergiefunktionals für dispersive Wechselwirkungen fällt auch hier das Helmholtzfunktional für dipolare Wechselwirkungen im homogenen Fall nicht auf die Formulierung innerhalb der PCP-SAFT Zustandsgleichung zurück. Daher kommt auch hier eine Korrektur zum Tragen, die mit der lokalen Dichteapproximation (LDA) ausgewertet wird. Dieses angepasste Funktional bzw. dessen Ableitung spiegelt sich in der Form δFDD[ρ̂(r,ω)] δρ̂(r,ω) = δF dd[ρ̂(r,ω)] δρ̂(r,ω) + ( µdd,EoS(r)− µdd(r,ω) ) . (4.54) mit dem hochgestellten Index ”DD”wider. Der hochgestellte Index ”dd” bezieht sich auf das orientierungsabhängige Funktional aus Gleichung (4.53). ”dd,EoS” bezieht sich auf den Dipol-Dipol-Beitrag in der Kernphase aus der PCP-SAFT Zustandsgleichung, in dem alle Orientierungsabhängigkeiten ausintegriert sind. Die Potentiale µdd,EoS(r) und µdd(r,ω) hängen neben der Temperatur von den Dichten ρ(r) bzw. ρ̂(r,ω) ab. 4.2.5.1 Dipol-Dipol-Wechselwirkungspotential Das im vorherigen Kapitel vorgestellte Helmholtzenergiefunktional für dipolare Wech- selwirkungen hängt maßgeblich vom Dipolwechselwirkungspotential uDD ab. Abbil- dung 4.1 zeigt die Winkelabhängigkeiten für zwei Teilchen und deren Dipole. Die vier Winkel, die die Orientierung der beiden Dipolmomente charakterisieren, sind in zwei Vektoren ω = (θ, φ)t und ω ′ = (θ′ , φ′)t zusammengefasst. Die entsprechen- den dreidimensionalen Einheitsvektoren m(ω) = (sin θ cosφ, sin θ sinφ, cos θ)t und m(ω′) = (sin θ′ cosφ′, sin θ′ sinφ′, cos θ′)t beschreiben die Ausrichtung der Dipolmo- mente in Zylinderkoordinaten. Der Vektor r = (r cos ζ, r sin ζ, z)t bezeichnet den Rich- tungsvektor zwischen beiden Teilchen, ζ ist der Winkel in der x, y-Ebene in Zylinderko- 28 KAPITEL 4. DICHTEFUNKTIONALTHEORIE Abbildung 4.1: Winkeldefinitionen der Dipolmomente zweier Teilchen in Zylinderko- ordinaten. ordinaten. Das Wechselwirkungspotential zwischen beiden dipolaren Teilchen lautet uDD(r,ω,ω′) = µ2 r̂3 ( m(ω)m(ω′)− 3m(ω)rm(ω′)r r̂2 ) , (4.55) mit r̂ = |r| als Abstand zwischen beiden Teilchen und µ als dimensionsloses Dipolmo- ment. Nach Auswertung der Vektorprodukte ergibt sich uDD(r,ω,ω′) = µ2 r̂3 ( cos θ cos θ′ + cosφ cosφ′ sin θ sin θ′ + sin θ sin θ′ sinφ sinφ′ − 3 r̂2 ( z cos θ + r cos(ζ − φ) sin θ )( z cos θ′ + r cos(ζ − φ′) sin θ′ )) . (4.56) 29 5 Grundlagen der Nichtgleichgewichts- Thermodynamik an Grenzflächen 5.1 Entropieproduktion aus Wärme- und Stofftransport durch Grenzflächen Die Entropieproduktion σ stellt in der Thermodynamik den Anteil der Entropieände- rung eines Prozesses dar, der durch irreversible Abläufe bedingt ist. Betrachtet man mo- lekulare Prozesse durch eine planare Dampf-Flüssig-Grenzfläche mit der z-Koordinate senkrecht zur Grenzfläche, kann die Entropieproduktion der gesamten Grenzfläche (pro Flächeneinheit) beschreiben werden durch σs = ∫ zs,l zs,g σ(z)dz . (5.1) Hierbei ist σ(z) die lokale Entropieproduktion. Die Dicke der Grenzfläche ist definiert als ∆zs = zs,l − zs,g (5.2) mit zs,l als Anfangspunkt der Grenzfläche auf der Flüssigseite und zs,g entsprechend als Endpunkt auf der Gasseite. Eine gängige Formulierung der lokalen Entropieproduktion ist die Summe der Produkte aller Stromdichten Ji mit ihren konjugierten Triebkräften Xi für n unabhängige Stromdichten σ(z) = n∑ i=1 JiXi (5.3) Vernachlässigt man chemische Reaktionen innerhalb und den Transport elektrischer La- dung durch die Grenzfläche hindurch, lässt sich die lokale Entropieproduktion nach 30 KAPITEL 5. GRUNDLAGEN DER NICHTGLEICHGEWICHTS-THERMODYNAMIK AN GRENZFLÄCHEN Glavatskiy und Bedeaux49 zweckmäßig schreiben als σ(z) = n∑ i=1 Jξi · ( −∂(µ̃i/T ) ∂z ) + Je · ( − 1 T 2 ∂T ∂z ) . (5.4) Hierbei sind Jξi = ρi · vi die Stromdichte der i -ten Komponente, Je die Gesamtener- giestromdichte und µ̃i das chemische Potential der Komponente i. Die Gesamtenergie- stromdichte Je = J ′ q + n∑ i=1 Jξih̃i (5.5) beinhaltet die messbare Wärmestromdichte J ′ q und die partielle molare Enthalpie h̃i, die alle Stoffstromdichten Jξi in sich tragen. Der Superindex ′ wird hier aus Konsistenz- gründen zur Formulierung von Glavatskiy49 mitgeführt. Die Gesamtenergiestromdichte Je und alle Größen, die mit einer Tilde versehen sind, beinhalten den Anteil der kine- tischen Energie v2 2 . Die Verwendung der Gesamtenergiestromdichte Je ist von Vorteil, da diese über die Grenzfläche hinweg bei stationären Nicht-Gleichgewichtsbedingungen konstant ist. Als Konstitutivbeziehung, die den Zusammenhang zwischen Stromdichten und Trieb- kräften darstellt, wird ein linearer Ansatz in der Widerstandsform gewählt, gemäß Xk = m∑ j=1 RkjJj . (5.6) Hierbei steht Rkj für die Widerstandskoeffizienten der Grenzfläche gegen Stoff- und Wär- metransport. Gleichung (5.6) ist eine generische Gleichung, die die Stromdichten J mit ihren Triebkräften X in Verbindung setzen. Im vorliegenden Fall weist Jj auf Elemente des Vektors {Jξi , Je}, Xk beinhaltet die Elemente {( −∂(µ̃i/T ) ∂z ) , ( − 1 T 2 ∂T ∂z )} . Die Betrachtung wird fortan auf ein stationäres System reduziert mit zeitlich unverän- derlichen Stromdichten durch die Grenzfläche hindurch. Die Kraft-Fluss-Beziehung aus Gleichung (5.6) kann damit über die Grenzfläche integriert werden49 mit dem Ergebnis ( 1 T l − 1 T g ) = Re qqJe + n∑ i=1 Re qiJξi − ( µ̃l j T l − µ̃g j T g ) = Re jqJe + n∑ i=1 Re jiJξi . (5.7) Rqq definiert den Widerstand gegen Wärmetransport, Rij den Widerstand gegen Stoff- transport und Riq den gekoppelten Widerstand. Der hochgestellte Index e gibt die Zu- ordnung zur Gesamtenergiestromdichte Je an. In ingenieursbezogenen Anwendungen wird für gewöhnlich mit dem messbaren Wärme- KAPITEL 5. GRUNDLAGEN DER NICHTGLEICHGEWICHTS-THERMODYNAMIK AN GRENZFLÄCHEN 31 strom J ′ q gearbeitet. Kombiniert man die Gleichungen (5.7) und (5.5), so ergeben sich für die Flüssigphase ( 1 T l − 1 T g ) = R′,l qqJ ′,l q + n∑ i=1 R′,l qiJξi − [( µ̃l j T l − µ̃g j T g ) − h̃l j ( 1 T l − 1 T g )] = R′,l jqJ ′,l q + n∑ i=1 R′,l jiJξi (5.8) und für die Gasphase ( 1 T l − 1 T g ) = R′,g qqJ ′,g q + n∑ i=1 R′,g qi Jξi − [( µ̃l j T l − µ̃g j T g ) − h̃g j ( 1 T l − 1 T g )] = R′,g jqJ ′,g q + n∑ i=1 R′,g ji Jξi . (5.9) Die hochgestellten Indizes (′, l) und (′, g) weisen auf den Zusammenhang zur messba- ren Wärmestromdichte J ′ q auf der Flüssig- und Gasphase hin. Nach Onsager50 sind die Kreuzeinträge der Widerstandsmatrix symmetrisch, wodurch gilt Rqi = Riq Rij = Rji . (5.10) Setzt man die Koeffizienten der Gleichung (5.7) mit (5.8) und (5.9) gleich, so ergeben sich die folgenden Beziehungen zwischen den Transportwiderständen für die Flüssig- und Gasphase R′,g qq = Re qq (5.11) R′,g qi = Re qi + h̃g iR e qq (5.12) R′,g ji = Re ji + h̃g iR e jq + h̃g jR e qi + h̃g i h̃ g jR e qq (5.13) R′,l qq = Re qq (5.14) R′,l qi = Re qi + h̃l iR e gq (5.15) R′,l ji = Re ji + h̃l iR e jq + h̃l jR e qi + h̃l ih̃ l jR e qq . (5.16) Die Formulierung ”absoluter” Enthalpien wie in den vorangegangenen Gleichungen für die Widerstandskoeffizienten ist nicht ohne weiteres möglich. ”absolut” bezieht sich hier- bei nicht auf den tatsächlichen Absolutwert, der die inneren Energieanteile der Atomker- ne etc. enthält, sondern nur der absolute Wert bezüglich den Standardbildungswerten. Hierzu wäre die Angabe der spezifischen Standardbildungenthalpien erforderlich. Im 32 KAPITEL 5. GRUNDLAGEN DER NICHTGLEICHGEWICHTS-THERMODYNAMIK AN GRENZFLÄCHEN Hinblick auf die Verwendung des DFT-Formalismus in Kombination mit der PC-SAFT Zustandsgleichung ist ein Ausdruck zu wählen, der residuelle Größen beinhaltet. Daher wird die partielle molare Enthalpie als residuelle Enthalpie ohne den idealen Gasanteil angesetzt. Dies steht nicht im Widerspruch zur Herleitung der Integralbeziehungen zur Berechnung der Widerstandskoeffizieten, wie im folgenden Abschnitt näher diskutiert wird. 5.2 Integralbeziehungen zur Beschreibung von Transportwiderstandskoeffizienten an Grenzflächen Um für die im vorangegangenen Abschnitt dargestellten Widerstandskeoffizienten eine Berechnungsvorschrift zu erhalten, ist die Formulierung des lokalen Widerstandsprofils notwendig. Johannessen und Bedeaux51 definierten einen Transportwiderstandskoeffizi- enten ganz allgemein wie folgt R = ∫ zs,l zs,g r(z) dz . (5.17) Hierbei ist r(z) der lokale Transporwiderstandskoeffizient innerhalb der Grenzfläche. Glavatskiy und Bedeaux49 entwickelten den Zusammenhang für ein Multikomponen- tensystem zwischen den lokalen Transportwiderstandskoeffizienten im Bezug auf den messbaren Wärmestrom einerseits und auf die Gesamtenergiestromdichte andererseits reqq = r′qq (5.18) reqi = −r′qqh̃i − n−1∑ k=1 r′qkxk + r′qi , i = (1, ...n− 1) (5.19) reqn = −r′qqh̃i − n−1∑ k=1 r′qkxk (5.20) reji = r′qqh̃jh̃i + n−1∑ k=1 xk(r ′ kqh̃i + r′qkh̃j)− (r′jqh̃i + r′qih̃j) + n−1∑ k=1 n−1∑ l=1 r′klxkxl − n−1∑ k=1 xk(r ′ ki + r′jk) + r′ji , j, i = (1, ...n− 1) (5.21) KAPITEL 5. GRUNDLAGEN DER NICHTGLEICHGEWICHTS-THERMODYNAMIK AN GRENZFLÄCHEN 33 reni = r′qqh̃nh̃i + n−1∑ k=1 xk(r ′ kqh̃i + r′qkh̃n)− r′qih̃n + n−1∑ k=1 n−1∑ l=1 r′klxkxl − n−1∑ k=1 r′kixk , i = (1, ...n− 1) (5.22) renn = r′qqh̃ 2 n + h̃n n−1∑ k=1 xk(r ′ kq + r′qk) + n−1∑ k=1 n−1∑ l=1 r′klxkxl . (5.23) Betrachtet man ein Reinstoffsystem und setzt diese Gleichungen über Gleichung (5.17) in die Beziehungen für die Transportwiderstandskoeffizienten (5.11) bis (5.16) ein, so erhält man Integralbeziehungen zur Berechnung des gasseitigen und flüssigseitigen Transport- widerstandskoeffizienten (vgl. Johannessen und Bedeaux51) R′,g qq =R′,l qq = ∫ zs,l zs,g r′qq(z) dz (5.24) R′,g q1 = ∫ zs,l zs,g r′qq(z) · ( h̃g − h̃(z) ) dz (5.25) R′,l q1 = ∫ zs,l zs,g r′qq(z) · ( h̃l − h̃(z) ) dz (5.26) R′,g 11 = ∫ zs,l zs,g r′qq(z) · ( h̃g − h̃(z) )2 dz (5.27) R′,l 11 = ∫ zs,l zs,g r′qq(z) · ( h̃l − h̃(z) )2 dz . (5.28) Die Gleichungen (5.24) bis (5.28) zeigen auf, dass die Widerstandskoeffizienten gegen Stofftransport R′,g 11 und R′,l 11 sowie die gekoppelten Widerstandskoeffizienten R′,g q1 und R′,l q1 sich in Gas- und Flüssigphase voneinander unterscheiden. Für die Berechnung aller Widerstandskoeffizienten werden die Enthalpien h̃l und h̃g in den Kernphasen, das Ent- halpieprofil h̃(z) sowie das lokale Widerstandsprofil r′qq(z) über die Grenzfläche hinweg benötigt. Erweitert man die Betrachtung auf ein binäres Gemisch, so erhöht sich die Anzahl der Widerstandskeoffizienten auf sechs für jede Seite der Grenzfläche. Für die Gasseite der Grenzfläche ergibt dies R′,g qq = ∫ zs,l zs,g r′qq(z) dz (5.29) R′,g q1 = ∫ zs,l zs,g [ r′qq(z) · ( h̃g 1 − h̃1(z) ) + r′q1(z) ·x2(z) ] dz (5.30) R′,g q2 = ∫ zs,l zs,g [ r′qq(z) · ( h̃g 2 − h̃2(z) ) − r′q1(z) ·x1(z) ] dz (5.31) 34 KAPITEL 5. GRUNDLAGEN DER NICHTGLEICHGEWICHTS-THERMODYNAMIK AN GRENZFLÄCHEN R′,g 11 = ∫ zs,l zs,g [ r′qq(z) · ( h̃g 1 − h̃1(z) )2 + 2 · r′q1(z) ·x2(z) · ( h̃g 1 − h̃1(z) ) + r′11(z) · (x2(z)) 2 ] dz (5.32) R′,g 22 = ∫ zs,l zs,g [ r′qq(z) · ( h̃g 2 − h̃2(z) )2 − 2 · r′q1(z) ·x1(z) ( h̃g 2 − h̃2(z) ) + r′11(z) · (x1(z)) 2 ] dz (5.33) R′,g 12 = ∫ zs,l zs,g [ r′qq(z) · ( h̃g 1 − h̃1(z) )( h̃g 2 − h̃2(z) ) + r′q1(z) · [ x2(z) · ( h̃g 2 − h̃2(z) ) − x1(z) · ( h̃g 1 − h̃1(z) )] + r′11(z) ·x1(z) ·x2(z) ] dz . (5.34) Die entsprechenden Ausdrücke für die Flüssigseite sind R′,l qq = ∫ zs,l zs,g r′qq(z) dz (5.35) R′,l q1 = ∫ zs,l zs,g [ r′qq(z) · ( h̃l 1 − h̃1(z) ) + r′q1(z) ·x2(z) ] dz (5.36) R′,l q2 = ∫ zs,l zs,g [ r′qq(z) · ( h̃l 2 − h̃2(z) ) − r′q1(z) ·x1(z) ] dz (5.37) R′,l 11 = ∫ zs,l zs,g [ r′qq(z) · ( h̃l 1 − h̃1(z) )2 + 2 · r′q1(z) ·x2(z) · ( h̃l 1 − h̃1(z) ) + r′11(z) · (x2(z)) 2 ] dz (5.38) R′,l 22 = ∫ zs,l zs,g [ r′qq(z) · ( h̃l 2 − h̃2(z) )2 − 2 · r′q1(z) ·x1(z) ( h̃l 2 − h̃2(z) ) + r′11(z) · (x1(z)) 2 ] dz (5.39) R′,l 12 = ∫ zs,l zs,g [ r′qq(z) · ( h̃l 1 − h̃1(z) )( h̃l 2 − h̃2(z) ) + r′q1(z) · [ x2(z) · ( h̃l 2 − h̃2(z) ) − x1(z) · ( h̃l 1 − h̃1(z) )] + r′11(z) ·x1(z) ·x2(z) ] dz . (5.40) Eine wichtige Eigenschaft des von Bedeaux und Mitarbeitern49,52 erarbeiteten Formalis- mus ist, dass alle Widerstandskoeffizienten einzig durch drei lokale Widerstandsprofile r′qq(z), r ′ q1(z) und r′11(z) eindeutig festgelegt sind. Zur Berechnung aller Widerstandsko- effizienten sind weiterhin nur Größen notwendig, die aus der DFT berechenbar sind. So ist die Kenntnis der partiellen molaren Enthalpien h̃1(z) und h̃2(z) beider Phasen im KAPITEL 5. GRUNDLAGEN DER NICHTGLEICHGEWICHTS-THERMODYNAMIK AN GRENZFLÄCHEN 35 Gleichgewicht sowie die Konzentrationsprofile x1(z) und x2(z) erforderlich. Für die Beschreibung der notwendigen lokalen Widerstandsprofile haben Johannessen und Bedeaux53,54,55 folgenden Ansatz ähnlich eines Square-Gradient-Ansatzes vorge- schlagen r(z) = rg + ( rl − rg ) ρ(z)− ρg ρl − ρg + α ( ρl ρ(z) )β ( dρ(z) dz )2 . (5.41) Die ersten beiden Terme hiervon beschreiben den Übergang zwischen den lokalen Widerstandswerten in den Kernphasen gewichtet mit dem Verlauf der lokalen Dichte. Der dritte Term beschreibt den Widerstand an der Grenzfläche proportional zum quadrierten Gradienten der Dichte. Hierin kommen zwei anpassbare Modellparameter α und β zur Anwendung. Der Term ( ρl ρ(z) )β definiert die Position des Maximums im lokalen Widerstandsprofil. Für positive Werte für β geht dieser Term in der Flüssigphase gegen Eins. In der Gasphase wird der Term aufgrund des Unterschieds der Dichte in beiden Phasen groß. Der zweite Parameter α quantifiziert die Gesamthöhe des Maximums im Profil. Beide Formfaktoren sind im Allgemeinen als temperatur- und druckabhängig anzusehen. Jedoch zeigten Untersuchungen von Reinstoffen56, dass beide Faktoren als konstant im untersuchten Temperatur- und Druckbereich angesehen werden können. Daher sind in den in dieser Arbeit gezeigten Ergebnissen ebenfalls beide Modellparameter als Konstanten für das jeweilige Stoffsystem angenommen. Bei der Auswertung aller hergeleiteten Integralbeziehungen ist schließlich die Wahl der Integrationsgrenzen, also die Definition der Grenzflächendicke entscheidend. Pu- blikationen im Bereich der NEMD Simulationen57,58,59,60 schlagen ein Kriterium vor, um die Grenzflächendicke zu definieren. Dieses Kriterium kommt auch bei den hier durchgeführten Berechnungen zum Einsatz und wird in Abschnitt 6.4 näher beschrieben. 36 6 Berechnungsablauf und numerische Umsetzung 6.1 Grenzflächenspannung dipolarer Reinstoffsysteme Wegen der Vielzahl an Orts- und Winkelintegrationen, die im allgemeinen dreidimen- sionalen Fall ausgewertet werden müssen, wird hier nur der Fall einer eindimensionalen, planaren Grenzfläche betrachtet. Damit ist das System invariant im Bezug auf Trans- lation und Rotation in der x,y-Ebene. Die in Kapitel 4.2.5 dargestellten Ortsvektoren reduzieren sich zu r = (z, 0, 0)t und r′ = (z + z′, r′, ζ ′)t 27. Abbildung 6.1 zeigt die rele- vanten Koordinaten im reduzierten System. Die Intergration über den radialen Abstand r′ wird durch die Integration über den Gesamtabstand r̂ = |r− r′| zwischen zweier Teil- chen ersetzt. Sowohl die Dichte als auch die Orientierungsverteilungsfunktion hängen Abbildung 6.1: Zwei dipolare Teilchen in Zylinderkoordinaten. Die Komplexität zum allgemeinen, dreidimensioalen Fall ist stark reduziert. somit nur noch von z and θ ab. Die Ortsdichtefunktion erhält man durch Integration KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG 37 von ρ̂(z, θ) über θ ρ(z) = ∫ π 0 ρ̂(z, θ) sin θ dθ . (6.1) Die Orientierungsverteilungsfunktion ergibt sich aus dem Quotienten beider Dichten α(z, θ) = ρ(z, θ) ρ(z) . (6.2) Für isotrope Fluide reduziert sich die Orientierungsverteilungsfunktion zu einer Kon- stanten αiso(z, θ) = ∫ 2π 0 αiso(r,ω) dφ = 1 2 . (6.3) Die Dichte ergibt sich zu ρ̂iso(z, θ) = 1 2 · ρ(z) . (6.4) Mit diesen Vereinfachungen ergeben sich die Grundgleichungen der Dichtefunktional- thoerie aus Kapitel 4.2.5 neu. Das Großkanonische Potential ist Ω = F [ρ̂(z, θ)]− ∫ ∫ µ∗(z, θ)ρ̂(z, θ) sin θdθdz . (6.5) Im thermodynamischen Gleichgewicht gilt δΩ δρ̂(z,θ) ∣ ∣ ∣ eq = 0 und somit δF [ρ̂(z, θ)] δρ̂(z, θ) = µ∗(z, θ) . (6.6) Der ideale Gasanteil, dessen Funktionalableitung sowie die Helmholtzenergiedichte re- duzieren sich zu F ig[ρ̂(z, θ)] kT = x ρ̂(z, θ) · [ ln ( Λ̂3ρ̂(z, θ) ) − 1 ] sin θdθdz (6.7) δF ig[ρ̂(z, θ)]/kT δρ̂(z, θ) = ln ( Λ̂3ρ̂(z, θ) ) (6.8) f ig[ρ̂(z, θ)] kT =ρ̂(z, θ) · [ ln ( Λ̂3ρ̂(z, θ) ) − 1 ] , (6.9) mit Λ̂3 = Λ3 · 2. Der Helmholtzenergiebeitrag durch dipolare Wechselwirkungen lässt sich in den gegebe- 38 KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG nen Variablen komplett formulieren als F dd[ρ̂(z, θ)] = 1 2 ∫ ∞ −∞ ∫ ∞ −∞ ∫ ∞ |z′| ∫ π 0 ∫ π 0 ρ̂(z, θ)ρ̂(z + z′, θ′)ghs(r̂, ρ̄)r̂ (6.10) × [∫ 2π 0 ∫ 2π 0 ∫ 2π 0 UDD(r̂, z ′, ζ ′, θ, φ, θ′, φ′)dφdφ′dζ ′ ] ︸ ︷︷ ︸ UDD(r̂,z′,θ,θ′) dθdθ′dr̂dzdz′ . Die Abhängigkeiten des Dipol-Dipol-Wechselwirkungspotentials reduzieren sich auf UDD(r̂, z ′, ζ ′, θ, φ, θ′, φ′). Die Auswertung durch ein Computeralgebra Programm ergibt einen vereinfachten Ausdruck von uDD mit uDD(r̂, z ′, θ, θ′) = µ2 · 8π3 r̂3 ( cos θ cos θ′ − 3z2 r̂2 cos θ cos θ′ ) . (6.11) Der quadrierte Ausdruck der Gleichung (4.56) liefert u2 DD(r̂, z ′, θ, θ′) = 1 2r̂10 ·µ4π3 · ( 9 · (r̂4 + z′4) ) − cos(2θ′) · ( r̂4 + 12r̂2z′2 − 27z′4 ) + cos(2θ) · ( −r̂4 − 12r̂2z′2 + 27z′4 ) + 9 cos(2θ′) · ( r̂4 − 8r̂z′2 + 9z′4 ) , (6.12) welcher in Gleichung (4.52) eingesetzt wird. Auf die Auswertung des kubischen Aus- drucks wird verzichtet. Die Integrationen über φ, φ′ und ζ ′ können analytisch aus- geführt werden, so dass sich die Abhängigkeiten des Wechselwirkungspotentials auf UDD(r̂, z ′, θ, θ′) reduzieren lassen. Damit ergibt sich für die Funktionalableitung sowie die Helmholtzenergiedichte im ein- dimensionalen Fall δF dd[ρ̂(z, θ)] δρ̂(z, θ) = ∫ ∞ −∞ ∫ ∞ |z′| ∫ π 0 ρ̂(z + z′, θ′) ( ghs(r̂, ρ̄) + 1 2 ρ(z) ∂ghs(r̂, ρ̄) ∂ρ̄(r̂) ) r̂ × UDD(r̂, z ′, θ, θ′) sin θ′dθ′dr̂dz′ (6.13) fdd[ρ̂(z, θ)] = 1 2 ρ̂(z, θ) ∫ ∞ −∞ ∫ ∞ |z′| ∫ π 0 ∫ π 0 ρ̂(z + z′, θ′) ghs(r̂, ρ̄) r̂ × UDD(r̂, z ′, θ, θ′) sin θ′dθ′dr̂dz′ . (6.14) Um eine Iterationsvorschrift für die Dichte ρ̂(z, θ) zu formulieren, werden die beiden Sei- ten von Gleichung (4.40) anhand Gleichung (4.10) aufgeteilt. Der Anteil res,ø beinhaltet KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG 39 den Hartkugel-Beitrag, den Kettenbeitrag sowie den dispersiven Beitrag. Dies liefert δF [ρ̂(z, θ)] δρ̂(z, θ) = kT ln ( Λ̂3ρ̂(z, θ) ) ︸ ︷︷ ︸ δF ig[ρ̂(z,θ)] δρ̂(z,θ) + δF res,ø[ρ(z)] δρ(z) + δFDD[ρ̂(z, θ)] δρ̂(z, θ) . (6.15) Wird nun die lokale Dichteapproximation aus Gleichung (4.54) eingesetzt, so ist ersicht- lich, dass sich in den Kernphasen die dipolaren Funktionalanteile aufheben. Im Bereich rund um die Grenzfläche überlagert der dipolare Funktionalanteil den Kernphasenbei- trag. Das chemische Potential ergibt sich durch Anwendung der lokalen Dichteapproximation laut Gleichung (4.54) zu µ(z, θ) = µig(z, θ) + µres,ø(z) + ✘ ✘ ✘ ✘✘µdd(z, θ) + ( µdd,EoS(z)− ✘ ✘ ✘ ✘✘µdd(z, θ) ) . (6.16) Hierbei hebt sich der Funktionalanteil der Dipol-Dipol-Wechselwirkungen auf und es bleibt nur noch der dipolare Beitrag der PCP-SAFT Zustandsgleichung übrig. Im thermodynamischen Gleichgewicht gilt das chemische Potential als konstant über die Grenzfläche hinweg. Daher kann das chemische Potential an jeder Position ausge- wertet werden. Da die PCP-SAFT Zustandsgleichung das chemische Potential in den Kernphasen liefert, wird hier der Wert aus der Flüssigphase (b,L) als Referenz gewählt µ(z, θ) = µEoS b,L = kT ln ( Λ̂3ρ̂b,L ) ︸ ︷︷ ︸ µig b,L +µres,ø b,L + µdd,EoS b,L ︸ ︷︷ ︸ µres,EoS b,L . (6.17) Die Kombination aus den Gleichungen (6.15) und (6.17) ergibt die Iterationsvorschrift für das Dichteprofil ρ̂(z, θ) =ρ̂b,L · exp ( µres,ø b,L /kT + µdd,EoS b,L /kT − δF res,ø[ρ(r)]/kT δρ(z) − δF dd[ρ̂(z, θ)]/kT δρ̂(z, θ) − ( µdd,EoS(z)/kT − µdd(z, θ)/kT ) ) =ρ̂b,L · exp ( µres,EoS b,L /kT − δF res,ø[ρ(r)]/kT δρ(z) − δFDD[ρ̂(z, θ)]/kT δρ̂(z, θ) ) . (6.18) Die Ausführung in Form einer Picard-Iteration als direktes Substitutionsverfahren lie- fert schließlich unter Berücksichtigung eines festgelegten Toleranzkriteriums ein entspre- chend ausiteriertes Profil der Dichte über die Grenzfläche in Abhängigkeit des Winkels θ. Für das Orstdichteprofil muss abschließend noch über θ integriert werden. Für die numerische Umsetzung von Gleichung (6.18) ist es zweckmäßig, das Dipolmo- ment, die Temperatur sowie die Ortskoordinaten in dimensionsloser Form zu verwenden. 40 KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG Die dimensionslosen Werte sind µ∗2 = µ2 εσ3 (6.19) T ∗ = T · k ε (6.20) z∗ = z σ r̂∗ = r̂ σ (6.21) mit der Boltzmannkonstanten k, dem Energieparameter ε und dem Segmentdurchmes- ser σ. Sollte der Segmentdurchmesser σ in [Å], der Energieparameter ε k in [K] und das Di- polmoment µ in [D] gewählt sein, muss Gleichung (6.19) folgendermaßen umformuliert werden µ∗2 = µ2kB εσ3 104 1, 3807 Å 3 K D2 . (6.22) 6.2 Grenzflächenspannung von binären Dampf-Flüssig-Systemen Der Berechnungsablauf zur Berechnung der Grenzflächenspannung binärer Dampf- Flüssig-Gleichgewichte und Reinstoffsysteme unterscheidet sich in nur wenigen Punkten. Gross27 beschreibt die grundlegenden Schritte des DFT-Formalsimus, welche im Kern hier ebenfalls Anwendung finden. Die Ortskoordinate normal zur gedachten Grenzflä- che wird in 1000 Stützstellen unterteilt mit einem äquidistanten Abstand zwischen den Stützpunkten von 0,25 Å. Die einzelnen Schritte des Berechnungsablaufs sind im Fol- genden beschrieben: 1. Vorab der Gleichgewichtsberechnung und Auswertung der Helmholtzenergiefunk- tionale wird die radiale Paarverteilungsfunktion für zwei Ketten der Komponente 1 ghc11 (r̂, η̂, m11), zwei Ketten der Komponente 2 ghc22 (r̂, η̂, m22), zwei Ketten der Pseu- dokomponente 1,2 ghc12 (r̂, η̂, m12) sowie deren partiellen Ableitungen ∂ghc(r̂,η̂,mij) ∂η̂ ausgewertet. Die Ergebnisse werden für Stützstellen, die den Variablenraum {r̂, η̂} abdecken, abgelegt. Die spätere Auswertung erfolgt über eine kubische Spline- Interpolation zwischen den berechneten Stützstellen. 2. Über die Auswertung der PCP-SAFT Zustandsgleichung wird zunächst das Pha- sengleichgewicht gelöst und die Kernphasengrößen Dichte und Zusammensetzung für beide Phasen sowie chemisches Potential und Druck berechnet. 3. Zur Festlegung des Startdichteprofils und der Berechnung des Kapillarwellenanteils wird die kritische Temperatur der Mischung Tmix c benötigt. Diese wird anhand des Algorithmus von Heidemann und Khalil61 berechnet. KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG 41 4. Für die Iteration des Dichteprofils beider Komponenten wird ein Startdichteprofil für beide Komponenten angesetzt ρ0i (z) = 1 2 ( ρli − ρgi ) tanh ( 2.7 z σi ( Tmix c − T Tmix c )) + 1 2 ( ρli + ρgi ) . (6.23) 5. Für die anschließende Iteration der Dichteprofile werden die Funktionalableitungen der einzelnen Anteile berechnet. Die Auswertung der Integrationen erfolgt über die Trapezregel. 6. Die Dichteprofile beider Komponenten werden anhand der Iterationsgleichung (4.15) bis zu einer vorgegebenen Abweichtoleranz χ = 10−6 iteriert. Aufgrund des direkten Substitutionsverfahrens und der starken Gradienten an der Grenz- fläche muss die Dichteaktualisierung in jedem Schritt stark gedämpft werden. Als Kompromiss zwischen Rechenzeit und Stabilität der Rechnung ist ein Dämpfungs- faktor von 0,005 für die meistern Berechnungen geeignet. Die Berechnung erfolgt über die beiden Gleichungen ρk+1, ungedämpft i (z) = ρli · exp ( µres,L i /kT − δF res[ρk](z)/kT δρi(z ) (6.24) ρk+1 i (z) = ρki (z) + 0, 005 · ( ρk+1, ungedämpft i (z)− ρki (z) ) (6.25) mit dem Iterationsschritt k. Das Abbruchkriterium der Iteration ist folgenderma- ßen definiert χ = 1 M M∑ m=1 n∑ i=1 ( ρk+1, ungedämpft i − ρki )2 (6.26) mit der Gesamtzahl an Stützstellen M und der Gesamtzahl an Komponenten n. 7. Die Berechnung der Grenzflächenspannung erfolgt schließlich über die Gleichungen (4.7) und (4.8). 6.3 Grenzflächenspannung von binären Flüssig-Flüssig-Systemen Für die Berechnung der Grenzflächenspannung von Systemen mit einer Flüssig-Flüssig- Entmischung ist das generelle Vorgehen gleich dem eines Dampf-Flüssig-Gleichgewichts. Bei der numerischen Umsetzung wurden lediglich wenige Parameter angepasst. Zum einen wurde der Diskretisierungsschritt in der Ortskoordinate auf 0,1 Å verrin- gert, um die starken Gradienten an der Grenzfläche besser abbilden zu können. Flüssig- Flüssig-Systeme weisen eine schmalere Grenzfläche im Vergleich zu Dampf-Flüssig- 42 KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG Systemen auf, wodurch die Dichtegradienten steiler sind (vgl. Ergebnisse in den Kapiteln 7.2 & 7.3). Des Weiteren wurde das Startprofil für die Iteration der Dichte geringfügig geändert. Die Temperaturabhängigkeit wurde durch eine Konstante ersetzt, da die Ska- lierung mit der kritischen Entmischungstemperatur kein geeignetes Konzept ist. Das angepasste Startdichteprofil für zwei Flüssigphasen L1 und L2 lautet ρ (0) i (z) = 1 2 ( ρL1i − ρL2i ) tanh ( 0.6 z σi ) + 1 2 ( ρL1i + ρL2i ) . (6.27) Bei der Berechnung der Grenzflächenspannung wurde auf die Kapillarwellenkorrektur (vgl. Gleichung (4.8)) verzichtet, da dieser Kapillarwellenansatz die universellen kriti- schen Skalierungsgesetze von Dampf-Flüssig Systemen extrapoliert, was sich nicht auf Flüssig-Flüssig Entmischungen übertragen lässt, weil dort andere Ordnungsparameter maßgeblich sind. 6.4 Transportwiderstände In diesem Abschnitt werden die Kernpunkte des Berechnungsablaufes aller Transport- widerstandskoeffizienten sowohl für Reinstoffsysteme als auch für binäre Mischungen erläutert. Annahmen und Einschränkungen sind an den jeweiligen Stellen angegeben, an denen sie zur Anwendung kommen. 6.4.1 Reinstoffsysteme Der Berechnungsablauf für ein Reinstoffsystem lässt sich in sieben Schritte untergliedern: 1. Mittels der PC-SAFT Zustandsgleichung wird das Phasengleichgewicht bestimmt und die notwendigen Kernphaseneigenschaften Dichte (ρg, ρl), chemisches Poten- tial (µg = µl) sowie molare Enthalpie (hg, hl) beider Phasen im Gleichgewichtszu- stand bei vorgegebener Temperatur und vorgegebenem Druck ermittelt. 2. Für diskretisierte Ortskoordinaten wird die Iterationsgleichung ρi(r) = ρLi · exp ( µres,L i /kT − δF res[ρk]/kT δρi(r) ) ∀ i (6.28) genutzt, um Dichteprofile ρ(z) zwischen den Kernphasen über den DFT- Formalismus zu berechnen (vgl. Abschnitt 6.2). KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG 43 3. Über die thermodynamische Beziehung h(z) = µ+ T · s(z) = µ+ T · ( ∂f̂(z) ∂T ) V,N (6.29) wird das molare Enthalpieprofil berechnet. Hierbei sind N die Molekülzahl, s(z) das auf Moleküle bezogene Entropieprofil, f̂(z) = f(z)/ρ(z) das molekulare Helm- holtzenergieprofil und f(z) die Helmholtzenergiedichte. Die partielle Ableitung bei konstantem Volumen und konstanter Molekülzahl in Gleichung (6.29) wird nume- risch über ein kubisches Polynom ausgewertet ( ∂f̂(z) ∂T ) V,n = (6.30) ( f̂(z, T − 2 ·∆T )− 8 · f̂(z, T −∆T ) + 8 · f̂(z, T +∆T ) + f̂(z, T + 2 ·∆T ) 12 ·∆T ) V,n . Die gewählte Temperaturdifferenz beträgt ∆T = 0, 5K. 4. Um die Grenzflächendicke abzuschätzen, wird das von Johannessen et al.56 defi- nierte Kriterium verwendet κ = ∣ ∣ ∣ ∣ ∣ ρi ( zs,l(κ) ) − ρli ρli ∣ ∣ ∣ ∣ ∣ κ = ∣ ∣ ∣ ∣ ρi (z s,g(κ))− ρgi ρgi ∣ ∣ ∣ ∣ (6.31) Hierbei steht κ für die Abweichung zwischen lokaler Dichte und Dichte in der Kernphase. Johannessen et al. verglichen für die Definition von κ verschiedene Be- rechnungen mit NEMD Simulationsdaten. Hieraus ergab sich ein geeigneter Wert von κ = 0, 01. Mit diesem Kriterium lassen sich Anfangs- und Endpunkt der Grenzfläche zs,l und zs,l bestimmen. 5. Neben dem Enthalpieprofil h(z) muss das lokale Widerstandsprofil gegen rei- nen Wärmetransport r′qq(z) berechnet werden. Für die Auswertung des Square- Gradient-Ansatzes r(z) = rg + ( rl − rg ) ρ(z)− ρg ρl − ρg + α ( ρl ρ(z) )β ( dρ(z) dz )2 (6.32) müssen zunächst über die thermische Leitfähigkeit λ die Größen r′,lqq und r′,gqq be- stimmt werden. Nach Groot et al.62 gilt r′qq = 1 λ ·T 2 . (6.33) 44 KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG Falls NEMD Simulationsdaten für das betrachtete System vorliegen, können bei- de lokalen Widerstandskoeffizienten über Termperaturprofile ausgewertet werden. Hierzu wird folgende Differentialgleichung63 r′qq(z) = − 1 T (z)2J ′ q(z) dT (z) dz , (6.34) verwendet. 6. Die anpassbaren Modellparameter α und β in Gleichung (5.41) werden an NEMD Simulationsdaten angepasst. 7. Um alle Widerstandskoeffizieten der Grenzfläche zu erhalten, werden die in Ab- schnitt 5.2 aufgezeigten Integralbeziehungen ausgewertet. 6.4.2 Multikomponentenmischungen Die allgemeine Herangehensweise zur Berechnung der Widerstandskoeffizienten binärer Mischungen ist ähnlich der für Reinstoffsysteme. In den folgenden Schritten wird auf die Abweichungen des Berechnungsablaufes näher eingegangen: 1. Für die Berechnung der Gleichgewichtsdaten Dichte ρi, chemisches Potential µi, partielle molare Enthalpie hi sowie Zusammensetzung xi in den Kernpha- sen bei gegebener Temperatur und Druck muss zunächst der binäre Interak- tionsparameter k12 der Mischung bestimmt werden. Die Anpassung erfolgt an experimentelle Dampf-Flüssig-Gleichgewichtsdaten. Über die Berthelot-Lorenz- Kombinationsregeln werden die Kreuzparameter für den Segmentdurchmesser σ12 und den Energieparameter ε12 berechnet. σ12 = 1 2 · (σ11 + σ22) ε12 = √ ε11 · ε22 · (1− k12) (6.35) 2. Die Berechnung der molekularen Dichteprofile ρ1(z) und ρ2(z) beider Komponen- ten erfolgt analog zu den Ausführungen in Abschnitt 6.2. 3. Für ein Multikomponentensystem gilt allgemein für die partielle molare Enthalpie einer Komponente i hi(z) = µi + T · si(z) = µi + T · ( ∂fi(z) ∂T ) V,nj 6=i (6.36) Hierbei sind si die partielle molare Entropie und fi die partielle molare Helmholtz- energie. Alle Größen werden hier als partielle molare Größen bezeichnet, obwohl es KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG 45 sich streng genommen um Größen handelt, die auf die Molekülzahl bezogen sind. Um die jeweilige partielle molare Größe berechnen zu können, wird auf die wesent- liche Eigenschaft einer partiellen molaren Größe z zurückgegriffen z = n∑ i=1 xizi . (6.37) xi ist der Stoffmengenanteil und zi die partielle molare Größe einer Komponente i. Sowohl der Kettenbeitrag als auch der Beitrag durch dispersive Wechselwirkungen der Helmholtzenergie sind bereits als Summe über alle Komponenten n formuliert, was die eindeutige Formulierung der jeweiligen partiellen molaren Größe daraus vereinfacht. Die exakten Ausdrücke sind in Anhang B aufgeführt. Die partielle Ableitung ( ∂fi(z) ∂T ) V,Nj 6=i der partiellen molaren Helmholtzenergie fi wird wie in Schritt 3 im vorangegangenen Abschnitt numerisch berechnet. 4. Das κ-Kriterium (Gleichung (6.31)) für die Grenzflächendicke wird für alle Kom- ponenten ausgewertet. Die Komponente mit der breiteren Grenzfläche definiert den Anfangs- und Endpunkt zs,l und zs,g für das System. 5. Zusätzlich zum lokalen thermischen Widerstandsprofil r′qq(z) werden das lokale Widerstandsprofil gegen reinen Stofftransport r′11(z), das Profil des lokalen, ge- koppelten Widerstands r′q1(z) sowie die Konzentrationsprofile xi(z) benötigt. Die Konzentrationsprofile lassen sich einfach über die Beziehung xi(z) = ρi(z) ∑2 j=1 ρj(z) (6.38) berechnen. Die Kernphasenwerte des Widerstandskoeffizienten gegen Stofftransport und des gekoppeltenWiderstandskoeffizienten lassen sich über die folgenden Beziehungen64 berechnen r′11 = R ( 1 x1 + ( ∂ lnϕ1 ∂x1 ) T,p ) Dρ2 (6.39) r′q1 = − STρ1x 2 2R ( 1 x1 + ( ∂ lnϕ1 ∂x1 ) T,p ) λρ2 . (6.40) Hierbei sind R die universelle Gaskonstante, ϕ1 der Fugazitätskoeffizient, D der Diffusionskoeffizient, ST der Soretkoeffizient und λ der thermische Leitkoeffizient von Komponente 1. Die partielle Ableitung ( ∂ lnϕ1 ∂x1 ) T,p wird über die Auswertung 46 KAPITEL 6. BERECHNUNGSABLAUF UND NUMERISCHE UMSETZUNG des Zentraldifferenzenquotienten berechnet. 6. Für jedes Profil des lokalen Transportwiderstandskoeffizienten müssen die Modell- parameter α und β unabhängig voneinander an entsprechende Widerstandskoeffi- zienten aus NEMD Simulationen angepasst werden. 7. Um alle Widerstandskoeffizienten der Grenzfläche zu erhalten, werden die in Ab- schnitt 5.2 aufgezeigten Integralbeziehungen ausgewertet. 47 7 Ergebnisse 7.1 Grenzflächeneigenschaften dipolarer Reinstoffe Bisherige Studien zur Berechnung von Grenzflächeneigenschaften realer, stark polarer Stoffe nutzen Funktionale, die ein isotropes Orientierungsverhalten voraussetzen. Hierbei zeigt sich eine systematische Abweichung zwischen berechneten und gemessenen Daten, welche auf die Vernachlässigung anisotroper Effekte zurückzuführen ist. Lediglich für Modellfluide wie das Stockmeyer Fluid mit punktförmigen Dipolmomenten wurden er- weiterte Funktionale zur Beschreibung anisotroper, dipolarer Effekte vorgestellt48,65,66,67. Das in dieser Arbeit betrachtete Helmholtzenergiefunktional soll durch die Berücksichti- gung der radialen Paarverteilungsfunktion reale, stark dipolare Fluide abbilden können. Daher wurde für Modellrechnungen der Stoff Aceton als Stellvertreter für kurzkettige, stark dipolare Moleküle betrachtet. Die Reinstoffparameter von Aceton für die PCP-SAFT Zustandsgleichung17 sind in Ta- belle 7.1 aufgelistet. Die Berechnungsbox wurde auf eine Länge von 40σ festgelegt, wobei Tabelle 7.1: PCP-SAFT Reinstoffparameter für Aceton17. mi σi / Å ǫi/k / K µi / Da 2,745 3,274 232,99 2,88 a D= 3, 3356 · 10−30C ·m ein Diskretisierungsschritt im Ort von ∆z = 0, 02 ·σ zur Anwendung kam. Die Winkelin- tegration erfolgte diskretisiert in 15◦ Abschnitten, um den Rechenaufwand vertretbar zu halten. Außerdem wurde ein Abschneideradius für die Auswertung des dipolaren und dispersiven Wechselwirkungspotentials von rc = 9σ gewählt. Um die Auswirkungen durch die Einführung des zusätzlichen Freiheitsgrads der Orien- tierung der Moleküle zu erkennen, wurden drei Berechnungen durchgeführt. In der ersten Berechnung wurde für das Fluid eine isotrope Orientierungsverteilung vorgegeben, mit α = αiso = 1 2 = konst. (7.1) 48 KAPITEL 7. ERGEBNISSE In der zweiten Berechnung wurde die Orientierungsverteilungsfunktion für 50 Iterati- onsschritte konstant gehalten, bevor diese für die restlichen Iterationsschritte ebenfalls freigegeben und in die Minimierung des Grosskanonischen Potentials Ω miteinbezogen wurde. In der dritten Berechnung wurden sowohl die ortsbezogene Dichte als auch die Orientierungsverteilungsfunktion von Anfang an simultan iteriert. Abbildung 7.1 zeigt die Entwicklung der nicht-konvergierten Gre