Kopplung des Einzelpartikel- und des Zwei-Kontinua-Verfahrens fur die Simulation von Gas-Feststoff-Stromungen [PDF]

Zitiervorschau

Kopplung des Einzelpartikel- und des Zwei-Kontinua-Verfahrens fu ¨ r die Simulation von Gas-Feststoff-Str¨ omungen

Zur Erlangung des akademischen Grades eines Dr.-Ing. vom Fachbereich Bio- und Chemieingenieurwesen der Universit¨at Dortmund genehmigte Dissertation vorgelegt von Dipl.-Ing. Christof Gru ¨ ner aus Tarnowitz Tag der m¨undlichen Pr¨ufung: 22. Oktober 2004 1. Gutachter: Prof. Dr. Karl Strauß 2. Gutachter: Prof. Dr. Herbert Koch Dortmund 2004

Berichte aus der Strömungstechnik

Christof Grüner

Kopplung des Einzelpartikel- und des Zwei-Kontinua-Verfahrens für die Simulation von Gas-Feststoff-Strömungen

.

D 290 (Diss. Universität Dortmund)

Shaker Verlag Aachen 2004

Bibliografische Information der Deutschen Bibliothek Die Deutsche Bibliothek verzeichnet diese Publikation in der Deutschen Nationalbibliografie; detaillierte bibliografische Daten sind im Internet über http://dnb.ddb.de abrufbar. Zugl.: Dortmund, Univ., Diss., 2004

.

Copyright Shaker Verlag 2004 Alle Rechte, auch das des auszugsweisen Nachdruckes, der auszugsweisen oder vollständigen Wiedergabe, der Speicherung in Datenverarbeitungsanlagen und der Übersetzung, vorbehalten. Printed in Germany.

ISBN 3-8322-3458-6 ISSN 0945-2230 Shaker Verlag GmbH • Postfach 101818 • 52018 Aachen Telefon: 02407 / 95 96 - 0 • Telefax: 02407 / 95 96 - 9 Internet: www.shaker.de • eMail: [email protected]

Danksagung Die vorliegende Arbeit entstand im Rahmen meiner T¨atigkeit als wissenschaftlicher Angestellter am Lehrstuhl Energieprozeßtechnik & Str¨omungsmechanik des Fachbereiches Biound Chemieingenieurwesen der Universit¨at Dortmund. Herrn Prof. Dr. K. Strauß danke ich herzlich f¨ ur die Unterst¨ utzung dieser Arbeit und f¨ ur das mir entgegengebrachte Vertrauen. Mit seinen Anregungen und Ratschl¨agen hat er maßgeblich zum Gelingen der Arbeit beigetragen. ¨ Danken m¨ochte ich auch Herrn Prof. Dr. H. Koch f¨ ur die freundliche Ubernahme des Korreferates. Frau Prof. Dr. G. Sadowski und Herrn Prof. Dr. P. Walzel danke ich f¨ ur die ¨ Mitwirkung an der Pr¨ ufung sowie Herrn Prof. Dr. U. K¨oster f¨ ur die Ubernahme des Kommissionsvorsitzes. Den studentischen Hilfskr¨aften, Studien- und Diplomarbeitern danke ich f¨ ur ihren vorbildlichen Einsatz bei der Bearbeitung der gestellten Themen. Den Mitarbeitern und Kollegen des Lehrstuhls Energieprozeßtechnik & Str¨omungsmechanik sei f¨ ur ihre Herzlichkeit und Hilfsbereitschaft gedankt. Mein ganz besonderer Dank geb¨ uhrt meinen Eltern, meinen Geschwistern und deren Familien f¨ ur den langj¨ahrigen R¨ uckhalt, das große Interesse, mit dem sie meinen Werdegang verfolgt haben, sowie die stets vorbehaltlose Unterst¨ utzung meiner Ziele. Meiner Frau Silvia danke ich f¨ ur ihre liebevolle und verst¨andnisvolle Unterst¨ utzung. Ihr sonniges Gem¨ ut und ihre Geduld waren mir stets eine Quelle meiner Bem¨ uhungen.

Meinen Eltern in Dankbarkeit

Zusammenfassung In der hier vorliegenden Arbeit wird ein gekoppeltes Verfahren f¨ ur die Simulation von Gas-Feststoff-Str¨omungen entwickelt. Als Basis f¨ ur die Kopplung dienen zwei etablierte Simulationsverfahren. Dies ist zum einen das deterministische Zwei-Kontinua-Verfahren, bei dem sowohl die Gasphase als auch der dispers verteilte Feststoff als koexistierende Kontinua aufgefasst werden. Zum anderen wird das stochastische Einzelpartikel-Verfahren eingesetzt, bei dem die Beschreibung der Bewegung der Feststoffteilchen mit Hilfe einer detaillierten Einzelpartikelbetrachtung vorgenommen wird. Die Aufteilung der Str¨omungsgeometrie in stochastische und deterministische Zonen erfolgt unter Verwendung eines Gebietszerlegungskriteriums, das zwei Bedingungen miteinander kombiniert, welche die G¨ ultigkeit des Zwei-Kontinua-Verfahrens pr¨ ufen. Zu diesem Zweck wird ein funktionaler Zusammenhang entwickelt, der den so genannten KTGFDeformationsanteil beschreibt. Mit dessen Hilfe ist die Beschreibung von Deformationen einer Verteilungsfunktion der Partikelgeschwindigkeiten m¨oglich, die aus den Gradienten der granularen Temperatur und der Partikelgeschwindigkeit resultieren. Bei diesen Deformationen handelt es sich um Abweichungen der Verteilungsfunktion von der Maxwell-Verteilung, die in der Herleitung der beschreibenden Gleichungen f¨ ur das Zwei-Kontinua-Verfahren zugrundegelegt sind. Die Formulierung des KTGF-Deformationsanteils ist so gew¨ahlt, dass ein direkter Vergleich mit einer Bewertungsfunktion DEq,p , die weitergehende Deformationen einer Verteilungsfunktion von der Gleichgewichtsverteilung erfasst, erfolgen kann. Eine weitere Bedingung ist die G¨ ultigkeit der Kontinuumsannahme, die ebenfalls eine notwendige Voraussetzung f¨ ur das Zwei-Kontinua-Verfahren ist. Die logische Verkn¨ upfung beider Bedingungen erm¨oglicht eine Zuordnung eines der beiden numerischen Berechnungsverfahren zu einem bestimmten Bereich der Str¨omungsgeometrie. F¨ ur die Konstruktion des gekoppelten Simulationsverfahrens wird das entwickelte Kriterium in einen adaptiven Gebietszerlegungsalgorithmus integriert. Ausgehend von der Definition einer Randzone, welche die stochastische Zone von der deterministischen Zone trennt, kann einer r¨aumlichen Ver¨anderung der Zonen im Verlaufe der Simulation Rechnung getragen werden. Innerhalb des gekoppelten Simulationsverfahrens spielt der Informationsaustausch an der Grenzfl¨ache zwischen den beiden numerischen Verfahren eine entscheidende Rolle. Um eine konservative Behandlung der Bilanzgr¨oßen an dieser Grenzfl¨ache zu gew¨ahr¨ leisten, werden entsprechende Ubergangsbedingungen vorgestellt. F¨ ur die Berechnung der Feststoffphase mit dem gekoppelten Simulationsverfahren wird eine Zwei-Schritt-L¨osungstrategie vorgeschlagen. Diese sieht im ersten Schritt eine Initialisierung vor, in der auf effiziente Weise eine N¨aherungsl¨osung erzeugt wird. Im Anschluss daran erfolgt eine erste Gebietsaufteilung. Im zweiten Schritt findet die gekoppelte Berechnung der Feststoffphase auf der Basis der adaptiven Gebietszerlegung statt. Das vorgestellte gekoppelte Simulationsverfahren wird anhand von Simulationen einer GasFeststoff-Str¨omung in ausgew¨ahlten Str¨omungskonfigurationen getestet. Dabei kann eine ¨ gute Ubereinstimmung zwischen den Ergebnissen des gekoppelten Verfahrens und denen des Einzelpartikel-Verfahrens, das als bestm¨ogliche Approximation zu Referenzzwecken herangezogen wird, festgestellt werden.

Inhaltsverzeichnis 1 Einleitung 1.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Ziele dieser Arbeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 2 3

2 Simulation von Gas-Feststoff-Str¨ omungen 2.1 Grundlagen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Mathematische Berechnungsverfahren . . . . . . . . . . . . . . . . . . . . .

5 5 8

3 Formulierung der Bilanzgleichungen 3.1 Beschreibung der Gasphase . . . . . . . . . . . . . . . . . 3.1.1 Bewegungsgleichungen . . . . . . . . . . . . . . . . 3.1.2 Modellierung der Turbulenz . . . . . . . . . . . . . 3.2 Beschreibung der Feststoffphase . . . . . . . . . . . . . . . 3.2.1 Zwei-Kontinua-Verfahren . . . . . . . . . . . . . . . 3.2.2 Einzelpartikel-Verfahren . . . . . . . . . . . . . . . 3.3 Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Randbedingungen f¨ ur das Zwei-Kontinua-Verfahren 3.3.2 Randbedingungen f¨ ur das Einzelpartikel-Verfahren 3.4 Numerisches Berechnungsverfahren . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

. . . . . . . . . .

13 13 13 14 18 19 24 32 32 35 36

4 Analyse des Verhaltens der Feststoffphase 4.1 Das granulardynamische Gleichgewicht . . . . . . . . . . 4.2 Erkennung granulardynamischer Gleichgewichtszust¨ande 4.3 Granulare Transportkoeffizienten . . . . . . . . . . . . . 4.4 Schlussfolgerung . . . . . . . . . . . . . . . . . . . . . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

39 39 40 43 44

5 Gebietszerlegungskriterium 5.1 G¨ ultigkeitsbereich der Kontinuumsannahme 5.2 KTGF-Deformationsanteil DKT GF . . . . . . 5.3 Signifikanz-Niveau-Beschr¨ankung . . . . . . 5.4 Gebietszerlegungskriterium . . . . . . . . . . 5.5 Schlussfolgerung . . . . . . . . . . . . . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

47 47 49 55 56 56

6 Konstruktion eines adaptiven Simulationsverfahrens 6.1 Konzept der adaptiven Verfahrenskopplung . . . . . . . 6.2 Informationsaustausch an den Zonengrenzen . . . . . . 6.2.1 Korrektur der Widerstandskraft . . . . . . . . . 6.2.2 Modellierung der Geschwindigkeitskovarianz . . ¨ 6.2.3 Formulierung der Ubergangsbedingungen . . . . 6.3 Beschreibung der L¨osungstrategie . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

59 59 62 62 66 68 74

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

i

Inhaltsverzeichnis 6.4 6.5

Modifikation des Mehrgitterverfahrens . . . . . . . . . . . . . . . . . . . . Schlussfolgerung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7 Analyse des gekoppelten Simulationsverfahrens 7.1 Simulation einer Gas-Feststoff-Str¨omung . . . . . . . . . . . . . . . . . . 7.1.1 Zylinderumstr¨omung . . . . . . . . . . . . . . . . . . . . . . . . . 7.1.2 Umstr¨omung eines einseitigen Hindernisses . . . . . . . . . . . . . 7.2 Untersuchung der L¨osungsqualit¨at des gekoppelten Simulationsverfahrens 7.2.1 Zylinderumstr¨omung . . . . . . . . . . . . . . . . . . . . . . . . . 7.2.2 Umstr¨omung eines einseitigen Hindernisses . . . . . . . . . . . . . 7.3 Schlussfolgerung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Zusammenfassung und Ausblick

ii

76 79

81 . 81 . 82 . 88 . 91 . 91 . 96 . 100 103

Variablenverzeichnis Lateinische Buchstaben DKT GF D DEq Ep Fr I Ip Kgp Kn L M Ma F N P R Re Rep ReS S St T T TR Wg Y+ c cd ckt cM cR dp e f fp fw g

KTGF-Deformationsanteil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Abmessung des Str¨omungsgebietes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m] Bewertungsfunktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] kinetische Energie eines Partikels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [J] Froude-Zahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Einheitstensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Tr¨agheitsmoment einer Kugel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [kg m2 ] Geschwindigkeitskovarianz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m2 /s2 ] Knudsen-Zahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] L¨angenmaßstab . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m] Momentenvektor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [Nm] Mach-Zahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Kraftvektor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [N] Partikelanzahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Wahrscheinlichkeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Radius . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m] Reynolds-Zahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Partikel-Reynolds-Zahl. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[–] Reynolds-Zahl der Scherstr¨omung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Kontrollvolumengrenze . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m2 ] Stokes-Zahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Spannungstensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [N/m2 ] thermodynamische Temperatur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[K] Tensor der Reynoldsspannungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [N/m2 ] mittlere Eintrittsgeschwindigkeit der Gasphase . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m/s] dimensionsloser Wandabstand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Fluktuationsgeschwindigkeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m/s] Widerstandsbeiwert. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[–] Korrekturfaktor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Auftriebsbeiwert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Rotationsbeiwert . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Partikeldurchmesser . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[m] Stoßzahl . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Verteilungsfunktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [verschieden] Reibungsbeiwert Partikel-Partikel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Reibungsbeiwert Partikel-Wand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Vektor der Erdbeschleunigung. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[m/s2 ] iii

Variablenverzeichnis g g0 k kgp
m, m mp n np p t u v x

Partikelgewicht . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] radiale Verteilungsfunktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] turbulente kinetische Energie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m2 /s2 ] Volumenkraft des Impulsaustausches zwischen Gas- und Partikelphase . . . [N/m3 ] freie Wegl¨ange . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m] statistische Momente . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [verschieden] Masse eines Partikels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [kg] Einheitsnormalenvektor einer Oberfl¨ache . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Anzahldichte der Partikelphase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [1/m3 ] Druck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [N/m2 ] Zeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [s] Vektor der Geschwindigkeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m/s] Vektor der Geschwindigkeit eines Einzelpartikels . . . . . . . . . . . . . . . . . . . . . . . . . . [m/s] Ortsvektor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m]

Griechische Buchstaben Γ Ω Φ β ε η γp κp µ µt νC ω φ ρ σ τ τC τK τp Θ ζ

Eintrittsfl¨ache. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[m2 ] Kontrollvolumen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m3 ] Korrekturpolynom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Widerstandsfunktion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [kg/m3 /s] Dissipationsrate der turbulenten kinetischen Energie. . . . . . . . . . . . . . . . . . . . .[m2 /s3 ] Kolmogorowsche L¨angenskala . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m] Dissipationsrate der granularen Temperatur . . . . . . . . . . . . . . . . . . . . . . . . . . [kg/m/s3 ] Leitf¨ahigkeit der Partikelphase. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[kg/m/s] Viskosit¨at . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [Pa s] turbulente Viskosit¨at . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .[Pa s] Kollisionsfrequenz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [1/s] Vektor der Rotationsgeschwindigkeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [1/s] Volumenanteil . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Dichte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [kg/m3 ] Standardabweichung . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [verschieden] Subzeitschritt . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [s] Kollisionszeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [s] Kolmogorov-Zeitskala . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [s] aerodynamische Partikelrelaxationszeit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [s] granulare Temperatur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [m2 /s2 ] Volumenviskosit¨at . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [Pa s]

Indizes Eq M W iv

bezeichnet den Gleichgewichtszustand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Maxwell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] an der Wand . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–]

Variablenverzeichnis g in i, j, k p ps rel

bezeichnet die Gasphase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] bezeichnet die Eintrittsparameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] ganzzahliger Laufindex . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] bezeichnet die Partikelphase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] simulierte Partikel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] relativ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–]

Hochgestellte Symbole (i)

Bezeichnung des Grades der Approximation zur Boltzmann-Gleichung. . . . . . . .[–]

Abk¨ urzungen DNS KTGF LES PSIC

Direkte Numerische Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Kinetische Theorie granularer Fluide . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Large Eddy Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–] Particle Source In Cell . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . [–]

v

Variablenverzeichnis

vi

1

Einleitung

Mehrphasenstr¨omungen treten in zahlreichen Anwendungen der industriellen Praxis auf. In der chemischen Verfahrenstechnik sind Mehrphasenstr¨omungen in vielen Grundoperationen wie beispielsweise Verdampfung, Kondensation, Mischen oder Trennen vertreten. Dabei lassen sich diese in vier Kategorien einteilen: Gas-Fl¨ ussig-, Gas-Feststoff-, Fl¨ ussig-Feststoffund Dreiphasenstr¨omungen. Die genannten Str¨omungskategorien k¨onnen in unterschiedlichen Konfigurationen auftreten. So stellt eine Blasens¨aule ein Gas-Fl¨ ussig-System dar, in dem die Fl¨ ussigkeit die kontinuierliche und die Gasblasen die disperse Phase darstellen. Ein Fl¨ ussigkeitszerst¨auber zielt auf die Dispersion feinster Fl¨ ussigkeitstropfen ab, wobei hier das Gas die kontinuierliche und die Tropfen die disperse Phase bilden. Als Beispiel f¨ ur eine weitere Konfiguration eines Gas-Fl¨ ussig-Systems kann eine geschichtete Rohrstr¨omung dienen. Hier bilden sowohl die Fl¨ ussigkeit als auch das Gas jeweils eine kontinuierliche Phase. Konkretisiert man die genannten Beispiele, so ist eine Unterscheidung in Ein- und Mehrkomponenten Mehrphasenstr¨omungen notwendig. Eine Dampf/Wasser Str¨omung stellt eine Gas-Fl¨ ussig-Str¨omung einer Komponente dar, w¨ahrend eine Luft/Wasser Str¨omung ein Beispiel einer mehrkomponentigen Gas-Fl¨ ussig-Str¨omung ist. Die Vielfalt und zugleich die Komplexit¨at der unterschiedlichen Str¨omungsformen stellen hohe Anspr¨ uche an den Ingenieur in Hinblick auf die Entwicklung, Auslegung und Optimierung zugeh¨origer Anlagenkomponenten. Lange Zeit stellten empirische Erfahrungen die Grundlage f¨ ur das Design derartiger verfahrenstechnischer Systeme dar. Das Fehlen einer leistungsf¨ahigen Messtechnik sowie der M¨oglichkeit einer Vorausberechnung wichtiger Einflußgr¨oßen machten Verfahrensverbesserungen außerordentlich schwierig. Erst die Fortschritte auf dem Gebiet der Messanalyse und die Bereitstellung hoher Rechnerleistungen haben dazu gef¨ uhrt, dass praxisrelevante Aufgabenstellungen deutlich effizienter bearbeitet werden k¨onnen. So erm¨oglicht die Entwicklung neuer Messtechniken eine zuverl¨assige Prozesskontrolle und die Ermittlung fundamentaler Prozessparameter. Wachsende Computerressourcen unterst¨ utzen den Einsatz moderner, numerischer Verfahren, die das Werkzeug eines Ingenieurs erg¨anzen. Insbesondere der anhaltende Fortschritt in der Entwicklung von Parallelrechnern, die aus mehreren tausend Prozessoren bestehen k¨onnen, erm¨oglicht Einblicke in bislang unerforschte Gebiete zu denen auch Aspekte der Mehrphasenstr¨omungen geh¨oren. Dabei liegen die Vorteile numerischer Simulationen in der schnellen und kosteng¨ unstigen Durchf¨ uhrung von Parameterstudien, wodurch zeitintensive und teure experimentelle Untersuchungen zumindest teilweise ersetzt werden k¨onnen. Zudem sind bestimmte Gr¨oßen einer experimentellen Untersuchung nicht zug¨anglich, so dass numerische Berechnungen die einzige M¨oglichkeit bieten, entsprechende Zusammenh¨ange zu untersuchen.

1

1 Einleitung

1.1

Motivation

Bei der theoretischen Beschreibung von zweiphasigen Str¨omungen haben sich zwei Verfahren etabliert, die sich in ihren Annahmen und Einsatzvoraussetzungen, vor allem aber in der Genauigkeit der L¨osung unterscheiden. W¨ahrend der Beschreibung der Gasphase ¨ jeweils gleiche Uberlegungen zugrunde liegen, ergeben sich die Unterschiede aus der Betrachtungsweise der dispersen Phase. Im so genannten Zwei-Kontinua-Verfahren wird die disperse Phase ausschließlich durch deren mittlere Eigenschaften beschrieben. Die L¨osung liefert keinerlei Informationen u ¨ber das Verhalten einzelner Elemente der dispersen Phase, insbesondere ist es nicht m¨oglich aus den zur Verf¨ ugung stehenden Informationen die Bahnkurve eines Partikels, Tropfens oder einer Gasblase in der Str¨omungsgeometrie anzugeben. Aufgrund der deterministischen Berechnung der zugrunde liegenden Transportgleichungen zeichnet sich das Verfahren durch eine hohe Effizienz im Hinblick auf die Rechengeschwindigkeit aus, liefert aber in Bereichen inhomogener Bewegung der dispersen Phase, beispielsweise vor Str¨omungshindernissen, unzutreffende Ergebnisse, die eine korrekte Vorausberechnung der Str¨omungsentwicklung verhindern. Einschr¨ankend ist weiterhin die Forderung, dass die Abst¨ande zwischen den Teilchen sehr viel kleiner sein m¨ ussen als die Abmessungen des Str¨omungsfeldes. Zudem erweist es sich als problematisch geeignete Randbedingungen f¨ ur einen weiten Anwendungsbereich zu formulieren. Dieses Problem tritt bei dem zweiten Verfahren nicht auf. Das so genannte EinzelpartikelVerfahren bei dem stochastische Methoden zum Einsatz kommen, liefert alle beschreibenden Zustandsgr¨oßen jedes einzelnen diskreten Partikels innerhalb der Geometrie. Als Partikel sind in diesem Zusammenhang sowohl Gasblasen, als auch Fl¨ ussigkeitstropfen und Feststoffteilchen zu verstehen. Wegen seiner detaillierten Modelle bietet dieses Verfahren die derzeit maximale Vorhersage-Genauigkeit f¨ ur die numerische Berechnung praxisrelevanter, zweiphasiger Str¨omungen. Die Formulierung der Randbedingungen erfolgt hier ebenfalls wesentlich detaillierter als bei dem Zwei-Kontinua-Verfahren und stellt keine Schwierigkeit dar. Dem steht als Nachteil der hohe Berechnungsaufwand gegen¨ uber, insbesondere dann, wenn Wechselwirkungen zwischen den Phasen ber¨ ucksichtigt werden m¨ ussen. In zahlreichen Anwendungen mit dispersen, zweiphasigen Systemen existieren Bereiche innerhalb des Str¨omungsgebietes, in denen das Verhalten der dispersen Phase durch das individuelle Verhalten jedes einzelnen Teilchens bestimmt wird. Da die Annahmen f¨ ur den Einsatz des Zwei-Kontinua-Verfahrens hier nicht eingehalten sind, ist die Anwendung des Einzelpartikel-Verfahrens zwingend erforderlich. In weiten Teilen der Str¨omungsgeometrie jedoch, wie das beispielsweise in Bereichen ungest¨orter Partikelbewegung der Fall ist, erlaubt das Verhalten der dispersen Phase den Einsatz des effizienten Zwei-KontinuaVerfahrens. Dabei kann in technisch relevanten Str¨omungskonfigurationen der Anteil der Geometrie, in dem das Zwei-Kontinua-Verfahren einsetzbar ist, deutlich u ¨ berwiegen. Hier wird ein Potential zur Rechenzeitersparnis erkennbar, wenn es gelingt die Str¨omungsgeometrie im Verlaufe der numerischen Berechnung in zwei Gebiete mit jeweils deterministischer und stochastischer Beschreibungsweise der dispersen Phase einzuteilen. Das rechenintensive Einzelpartikel-Verfahren kann dann auf den Bereich eingeschr¨ankt werden, in dem es unverzichtbar ist. Im verbleibenden Großteil des Rechengebietes kann auf das effiziente 2

1.2 Ziele dieser Arbeit Zwei-Kontinua-Verfahren zur¨ uckgegriffen werden. Somit kann die numerische Berechnung wesentlich zeitsparender durchgef¨ uhrt werden, ohne einen nennenswerten Verlust in der Genauigkeit in Kauf nehmen zu m¨ ussen.

1.2

Ziele dieser Arbeit

In der vorliegenden Arbeit soll f¨ ur Konfigurationen turbulenter Gas-Feststoff-Str¨omungen ein neuartiges, effizientes Verfahren zur Beschreibung der dispersen Phase entwickelt werden. Ziel ist es, das stochastische Einzelpartikel-Verfahren und das deterministische Zwei-Kontinua-Verfahren auf der Basis einer adaptiven Gebietszerlegung zu koppeln. Das Hauptaugenmerk liegt dabei auf der Entwicklung von Kriterien und Techniken welche die Konstruktion des kombinierten Verfahrens erm¨oglichen. Mit dessen Hilfe sollen die Vorteile des numerisch sehr effizienten Zwei-Kontinua-Verfahrens mit der hohen Genauigkeit des Einzelpartikel-Verfahrens kombiniert werden. Im folgenden Kapitel werden zun¨achst wichtige charakteristische Eigenschaften von ¨ Gas-Feststoff-Str¨omungen erl¨autert sowie ein Uberblick u ¨ber die numerischen Berechnungsmethoden zur Beschreibung von Zwei-Phasen-Str¨omungen gegeben. In Kapitel 3 werden die Bilanzgleichungen zur Beschreibung der Gas- und der Feststoffphase vorgestellt. Dabei wird der Formulierung der Bilanzgleichungen der Feststoffphase sowohl der kontinuierlicher Charakter des Zwei-Kontinua-Verfahrens als auch die diskrete Beschreibungsweise des Einzelpartikels-Verfahrens Rechnung getragen. Im Kapitel 4 wird auf Basis der kinetischen Theorie granularer Fluide das Verhalten der Feststoffphase analysiert. Die Entwicklung eines Gebietszerlegungskriteriums, dessen Einsatz eine Voraussetzung f¨ ur die Durchf¨ uhrbarkeit der Verfahrenskopplung darstellt, steht im Mittelpunkt von Kapitel 5. In dem darauf folgenden Kapitel 6 wird das Konzept f¨ ur die Entwicklung eines adaptiven Simulationsverfahrens dargelegt. Sowohl die notwendigen Techniken zur Umsetzung des neuen Verfahrens als auch die erforderliche Gesamt-L¨osungstrategie werden hier beschrieben und diskutiert. In Kapitel 7 schließlich wird die L¨osungsqualit¨at und die Rechenzeitersparnis des entwickelten Simulationsverfahrens analysiert. Dabei werden anhand von Simulationen in ausgew¨ahlten Str¨omungskonfigurationen die Ergebnisse der Berechnung des Bewegungsverhaltens der Feststoffphase, die mit dem Einzelpartikel- und dem gekoppelten Simulationsverfahren erhalten wurden, miteinander verglichen.

3

1 Einleitung

4

2

Simulation von Gas-Feststoff-Str¨ omungen

Die Vielzahl technischer Anwendungsbereiche f¨ ur Gas-Feststoff-Str¨omungen impliziert bereits, dass ein breiter Wertebereich f¨ ur die str¨omungstechnisch relevanten Parameter besteht. Dabei ist die Kenntnis u ¨ber die Auspr¨agung dieser Parameter von großer Bedeutung f¨ ur die theoretische Modellbildung. So nimmt beispielsweise die H¨ohe des Volumens, das von der Feststoffphase in einem Bezugsvolumen eingenommen wird, einen entscheidenden Einfluss auf die Auswahl der Modelle zur Beschreibung der Wechselwirkungsmechanismen innerhalb der dispersen Phase. In diesem Kapitel werden zun¨achst die wichtigsten Kennzah¨ len erl¨autert, die eine Gas-Feststoff-Str¨omung charakterisieren. Danach wird ein Uberblick u ¨ber die verschiedenen Verfahren zur Simulation von Gas-Feststoff-Str¨omungen gegeben sowie eine Eingrenzung des Bereiches der in dieser Arbeit untersuchten Konfigurationen vorgenommen.

2.1

Grundlagen

Gas-Feststoff-Str¨omungen besitzen eine Reihe charakteristischer Merkmale anhand derer eine Klassifizierung m¨oglich ist. Ein wichtiger Parameter ist der relative Volumenanteil φ, der den Anteil des Volumens angibt, den eine Phase der Suspension pro Einheitsvolumen einnimmt. Bezeichnet man den Volumenanteil der Feststoffphase mit φp und den der Gasphase mit φg , so gilt f¨ ur die Summe beider Anteile: φp + φg = 1.

(2.1)

Wird nun eine Phase als Kontinuum betrachtet, so sind nur die makroskopischen Eigenschaften derselben interessant. Sie werden durch eine Mittelung der lokalen mikroskopischen Ereignisse gebildet. Der Bereich auf den sich diese Mittelwertbildung erstreckt, wird Kontrollvolumen genannt. Der Gr¨oße des Kontrollvolumens der Mittelwertbildung kommt eine entscheidende Rolle zu. Um eine statistisch signifikante Mittelwertbildung sicherzustellen, muss das Volumen eine ausreichend große Anzahl von dispersen Elementen, hier zusammenfassend als Partikel bezeichnet, enthalten. So ist beispielsweise bei einer Anzahl von 104 Partikeln innerhalb eines Kontrollvolumens der Fehler, den man bei einer Mittelwertbildung macht, kleiner als 1% [Reif65]. F¨ ur ein Gas unter Normalbedingungen w¨are hierzu ein w¨ urfelf¨ormiges Volumen mit einer Kantenl¨ange von ca 1µm erforderlich [Cro82]. Da die Abmessungen der meisten Str¨omungskonfigurationen deutlich gr¨oßer sind als diese Skala, ist hier die Kontinuumsannahme f¨ ur die Gasmolek¨ ule gerechtfertigt. Betrachtet man jedoch die disperse Phase einer Gas-Feststoff-Str¨omung, so kann das erforderliche Kontrollvolumen aufgrund eines geringen Volumenanteils der Feststoffphase oder eines großen Eigenvolumens der Partikel schnell die charakteristischen Dimensionen der Geometrie erreichen, wenn der Fehler ebenfalls weniger als 1% betragen soll. In solchen F¨allen ist die 5

2 Simulation von Gas-Feststoff-Str¨omungen Kontinuumsannahme f¨ ur die Feststoffphase nicht mehr zutreffend. Daher ist im Einzelfall zu u ufen, ob eine Kontinuumsbetrachtung der dispersen Phase zul¨assig ist. ¨berpr¨ Ein weiterer wichtiger Parameter zur Klassifizierung von Gas-Feststoff-Str¨omungen ist die Stokes-Zahl St. Diese ist definiert als das Verh¨altnis aus der aerodynamischen Partikelrelaxationszeit τp und einer charakteristischen Zeitskala des Str¨omungssystems, die in turbulenten Gas-Feststoff-Str¨omungen h¨aufig aus dem so genannten Kolmogorov-Zeitmaß τK gebildet wird: St =

τp . τK

(2.2)

Die aerodynamische Partikelrelaxationszeit τp ist definiert als die Zeitspanne in der ein Partikel aus seiner Ruhelage heraus 63%, dies entspricht dem Wert 1 − exp(−1), der station¨aren Endgeschwindigkeit erreicht. Sie kennzeichnet damit das dynamische Antwortver¨ halten der Partikel bei Anderung der umgebenden Gasphasengeschwindigkeit. τp kann aus den Eigenschaften des Systems abgeleitet werden: τp =

ρp d2p . 18 µg

(2.3)

Dabei kennzeichnen ρp und dp die Dichte und den Durchmesser eines Partikels, w¨ahrend µg die molekulare Viskosit¨at der Gasphase darstellt. Das Kolmogorov-Zeitmaß τK gibt die Lebensdauer der kleinsten turbulenten Wirbel der Gasphase an. Sie stellt somit ein Maß f¨ ur das dynamische Verhalten der Gasphase dar. F¨ ur τK gilt:
µg . (2.4) τK = ρg ε Hier sind ρg die Dichte der Gasphase und ε die Dissipationsrate der turbulenten kinetischen Energie k, auf die in Kapitel 3 genauer eingegangen wird. W¨ahrend τp und τK die ¨ Zeitskalen f¨ ur die Anderung der Geschwindigkeit jeweils einer Phase angeben, beschreibt die Stokes-Zahl den Grad des Folgeverm¨ogens der Partikel in der Gasphase. F¨ ur St
1 ¨ erfolgt die Reaktion der Partikelgeschwindigkeit auf eine Anderung der Gasgeschwindigkeit nahezu ohne Zeitverzug. Die Bahnlinie eines Partikels stimmt dann in guter N¨aherung mit der Bahnlinie eines benachbarten Fluidelementes u ur St  ¨berein (vgl. Abbildung 2.1). F¨ ¨ 1 dagegen ist die Anderungsgeschwindigkeit der Gasphase so groß, dass die Partikel aufgrund ihrer Tr¨agheit kaum darauf reagieren k¨onnen. Das Folgeverm¨ogen der Partikel ist in diesem Fall gering, sie werden in ihrer Bewegung nur geringf¨ ugig von den Fluktuationen der Gasphase beeinflusst. Bei St ≈ 1 neigen die Partikel dazu den turbulenten Strukturen der Gasphase zu folgen. Aufgrund von Zentrifugaleffekten kommt es jedoch zur Kr¨ ummung ihrer Bahnlinien. Tang, der diese Thematik numerisch untersucht hat, beobachtete eine Aufkonzentrierung der Partikel an den R¨andern der Wirbelstrukturen [TWY92], die von Yang experimentell nachgewiesen werden konnte [YCC95]. Ein weiterer wichtiger Klassifizierungsparameter ist die Kollisionszeit τC . Diese kennzeichnet die mittlere Zeit, die vergeht, bis ein Partikel nach einer Kollision erneut mit einem 6

2.1 Grundlagen Bahnlinie eines Partikels

St
1

Bahnlinie der Gasphase

St1 St≈1

Abbildung 2.1: Einfluss der Stokes-Zahl auf die Partikelbewegung in turbulenten Strukturen der Gasphase. anderen Partikel kollidiert. τC kann aus der Kollisionsfrequenz νC berechnet werden, auf die in Kapitel 3 genauer eingegangen wird: τC =

1 . νC

(2.5)

Wird nun die aerodynamische Partikelrelaxationszeit τp auf die Kollisionszeit τC bezogen, so k¨onnen wichtige Aussagen u ¨ber die dominanten Phasenwechselwirkungen getroffen werden. Ist der Quotient τp
1 (2.6) τC so spricht man von einer durch aerodynamische Kr¨afte dominierten Gas-Feststoff-Str¨omung. Hier wird die Partikelbewegung maßgeblich von der Gasphase beeinflusst. Die Partikel ha¨ ben gen¨ ugend Zeit um auf Anderungen der Gasgeschwindigkeit zu reagieren, bevor eine Kollision stattfindet. In diesem Zusammenhang spricht man auch von einer verd¨ unnten Str¨omung [CST98]. Ist auf der anderen Seite τp  1, (2.7) τC so liegt eine kollisionsdominierte Str¨omung vor. Die Partikel kollidieren im zeitlichen Mit¨ tel mit einem anderen Partikel, bevor sie auf die Anderungsgeschwindigkeit der Gasphase reagieren k¨onnen. In diesem Parameterbereich verwendet man den Begriff einer dichten Gas-Feststoff-Str¨omung [CST98]. Zwischen der Kollisionszeit τC und dem Volumenanteil besteht ein einfacher Zusammenhang: steigt der Volumenanteil an, so verringert sich die Kollisionszeit τC aufgrund kleiner werdender Abst¨ande zwischen den Partikeln. Der Quotient aus τp und τC wird dementsprechend gr¨oßer. Nach einer Einteilung von Elghobashi ¨ findet bei einem Volumenanteil von φp = 10−3 der Ubergang von einer verd¨ unnten zur einer dichten Suspension statt [Elg94]. Bei geringeren Volumenanteilen d.h. beim Vorliegen einer verd¨ unnten Suspension k¨onnen die interpartikul¨aren Wechselwirkungen vernachl¨assigt werden. Eine Phasenwechselwirkung kann dagegen auch bei weit kleineren Volumenanteilen eine wichtige Rolle spielen. So ist in einem Volumenanteil-Bereich zwischen 7

2 Simulation von Gas-Feststoff-Str¨omungen φp = 10−6 und φp = 10−3 die R¨ uckwirkung der Partikel auf die Bewegung der Gasphase in die Modellierung einzubeziehen. In diesem Parameterbereich spricht man von einer Zwei-Wege-Kopplung. Dabei wird sowohl die Geschwindigkeit als auch die Intensit¨at der ¨ Gasphasenturbulenz von der Pr¨asenz der Feststoffphase beeinflusst. Die Anderung der Gasphasengeschwindigkeit resultiert aus der Impuls¨anderung der Feststoffphase. Die Intensit¨at der Gasphasenturbulenz kann in Anwesenheit der Partikel sowohl angefacht als auch ged¨ampft werden. So erh¨ohen Partikel in einem Gas-Feststoff-System mit St < 100 die Turbulenzdissipation, bei h¨oheren St-Zahlen dagegen wird die Turbulenzproduktion verst¨arkt [Elg94]. Bei Volumenanteilen φp < 10−6 liegt eine Ein-Wege-Kopplung vor, da lediglich die Gasphase einen Einfluss auf die Partikel aus¨ ubt. Erh¨oht man den Volumenanteil der dispersen Phase u ¨ber einen Wert von φp > 10−3 , liegt also eine dichte Suspension vor, so m¨ ussen auch Partikel-Partikel-Kollisionen in die Modellierung einbezogen werden. In diesem Zusammenhang spricht man von einer Vier-Wege-Kopplung. In der vorliegenden Arbeit werden sowohl verd¨ unnte als auch dichte Gas-FeststoffStr¨omungen untersucht. Die Einbeziehung von Str¨omungshindernissen in die Gesamtbetrachtung verursacht hohe lokale Unterschiede in der Auspr¨agung der einzelnen Str¨omungsparameter, wie beispielsweise des Volumenanteils. Deshalb m¨ ussen a priori alle Wechselwirkungsmechanismen zwischen der Gas- und der Feststoffphase in die Betrachtung einbezogen werden, um eine realit¨atsnahe Beschreibung derartiger Gas-Feststoff-Systeme zu gew¨ahrleisten.

2.2

Mathematische Berechnungsverfahren

Die Entwicklung mathematischer Berechnungsverfahren, die eine vollst¨andige Beschreibung von Gas-Feststoff-Str¨omungen erm¨oglichen, stellt seit langem eine Herausforderung an die Wissenschaft dar. Obwohl ein allgemein g¨ ultiger L¨osungsansatz bislang nicht verf¨ ugbar ist, existieren bereits mehrere Modellans¨atze, die eine Vorhersage des Bewegungsverhaltens einer Gas-Feststoff-Str¨omung erm¨oglichen und sich insbesondere in industriellen Anwendungen bew¨ahrt haben. Neben einer einphasigen Berechnung beider Phasen konnten sich vor allem die zweiphasigen Berechnungsverfahren etablieren. Zu den Letzten geh¨oren das Zwei-Kontinua-Verfahren und das Einzelpartikel-Verfahren. Beide Verfahren unterscheiden sich in der Beschreibung der Feststoffphase und werden zur Berechnung von verd¨ unnten sowie dichten Gas-Feststoff-Str¨omungen bis zu einem Volumenanteil von 10−1 eingesetzt. Bevor die Grundideen der zweiphasigen Berechnungsverfahren vorgestellt werden, soll kurz auf die Methode der einphasigen Beschreibung einer Zweiphasenstr¨omung eingegangen werden. Wie bereits im letzten Abschnitt erw¨ahnt, implizieren kleine St-Zahlen (St
1) ein sehr gutes Folgeverm¨ogen der Partikel gegen¨ uber der Gasphase. In einem solchen Fall k¨onnen beide Phasen als eine quasi Einphasenstr¨omung betrachtet werden. Hierzu werden die beschreibenden Gr¨oßen der Bewegung und der Stoffeigenschaften jeder Einzelphase zu Mischungsgr¨oßen zusammengefasst. Der Vorteil dieses Ansatzes besteht in der Reduktion der Anzahl der beschreibenden Differentialgleichungen des Systems. Aus diesem Grund besitzt das Verfahren eine hohe numerische Effizienz. Als nachteilig ist die mangelnde phy8

2.2 Mathematische Berechnungsverfahren sikalische Motivation bei der Definition einiger Mischungsgr¨oßen zu nennen, insbesondere dann, wenn das Folgeverm¨ogen der Partikel gegen¨ uber der Gasphase nachl¨asst. Desweiteren werden Phasenwechselwirkungen nicht erfasst, und Informationen u ¨ber das Bewegungsverhalten der dispersen Phase sind nicht verf¨ ugbar. Da f¨ ur das im Rahmen dieser Arbeit verfolgte Ziel eine detaillierte Beschreibung der dispersen Phase erforderlich ist, kann hier die einphasige Beschreibung keine Anwendung finden. Zwei-Kontinua-Verfahren Untersuchungsziel einer Gas-Feststoff-Str¨omung ist oftmals nicht die vollst¨andige, mikroskopische Beschreibung des gesamten Str¨omungsbereiches, sondern die Kenntnis der ortsund zeitaufgel¨osten mittleren, makroskopischen Gr¨oßen, wie der Geschwindigkeit, des Volumenanteils oder der Temperatur beider Phasen. Das Zwei-Kontinua-Verfahren, in der Literatur auch unter dem Namen Euler/Euler-Verfahren bekannt, behandelt Gas- und Feststoffphase als zwei koexistierende Kontinua und liefert so die mittleren Str¨omungseigenschaften beider Phasen. W¨ahrend die Kontinuumsannahme f¨ ur die Gasphase in den meisten F¨allen zutrifft - ausgenommen bei stark verd¨ unnten Gasen - sollte eine kontinuierliche Behandlung der Feststoffphase aufgrund des individuellen Bewegungsverhaltens der Partikel kritisch u uft werden. Deshalb ist die Wahl eines geeigneten Kontrollvolu¨berpr¨ mens f¨ ur die Kontinuumsbetrachtung eine wichtige Vor¨ uberlegung. Wie bereits im vorigen Abschnitt erw¨ahnt, sollte dieses zum einen eine ausreichende Anzahl von Partikeln enthalten, zum anderen sollte es klein genug gegen¨ uber den charakteristischen Abmessungen der Str¨omungsgeometrie sein. Um die mikroskopischen Effekte, die in einem Kontrollvolumen stattfinden, auf eine makroskopische Ebene zu u ¨bertragen, werden die zugrunde liegenden Transportgleichungen beider Phasen einem Mittelungsprozeß unterzogen, das in Kapitel 3 erl¨autert wird. Die Ber¨ ucksichtigung der Wechselwirkungen zwischen den Phasen erfolgt durch die Einf¨ uhrung lokaler, gemittelter Austauschterme. Die Hauptprobleme bei der Verwendung des Zwei-Kontinua-Verfahrens liegen einerseits in der Bestimmung der Transportkoeffizienten der Feststoffphase, wie beispielsweise der Partikelviskosit¨at, andererseits in der Formulierung geeigneter Randbedingungen f¨ ur die Feststoffphase. F¨ ur die Ermittlung der Transportkoeffizienten wurden in den letzten Jahrzehnten zahlreiche Ans¨atze entwickelt, die weitestgehend auf experimentelle Daten zur¨ uckgreifen und damit nur f¨ ur einen beschr¨ankten Parameterbereich g¨ ultig sind [CSL86, RE89]. Eine allgemein g¨ ultige M¨oglichkeit zur Berechnung der Transportkoeffizienten bieten analytische Ans¨atze, zu denen die weit verbreitete kinetische Theorie granularer Fluide (KTGF) geh¨ort. Zur Entwicklung dieser Theorie wurde die Analogie zwischen der Bewegung der Molek¨ ule eines Gases und der Bewegung der Partikel eines fluidisierten Feststoffs verwendet [CC70, Gid94]. Eine Einschr¨ankung bei der Verwendung dieses Ansatzes stellt die Annahme dar, dass die Partikelgeschwindigkeiten einer bestimmten Wahrscheinlichkeitsdichteverteilung unterliegen. Wird diese Verteilung durch ¨außere Einfl¨ usse, beispielsweise durch Hindernisse, stark deformiert, kann die KTGF nicht ohne Einschr¨ankungen verwendet werden [Kan03]. Eine weitere Schwierigkeit bei der Verwendung des Zwei-Kontinua-Verfahrens stellt die 9

2 Simulation von Gas-Feststoff-Str¨omungen Formulierung geeigneter Randbedingungen f¨ ur die Partikelphase dar. Wechselwirkungen zwischen der Feststoffphase und der Bewandung, wie beispielsweise inelastische Kollisionen oder Wandrauhigkeitseffekte, erschweren die Beschreibung solcher Vorg¨ange auf der makroskopischen Ebene. Um diesem Problem zu begegnen, wurden zahlreiche Modelle entwickelt, von denen jedoch viele problemspezifisch sind und daher keine Allgemeing¨ ultigkeit besitzen [Ind01]. Eine weitere Einschr¨ankung ergibt sich, wenn Partikelkollektive mit einer Gr¨oßenverteilung analysiert werden sollen. Die Einbeziehung einer polydispersen Partikelphase in die Simulation ist bei der Anwendung des Zwei-Kontinua-Verfahrens zwar m¨oglich, erfordert aber eine Einteilung des Partikelkollektivs in Gr¨oßenklassen [Boe96, Sya87]. Da jede Klasse als eine disperse Phase betrachtet werden muss, ist diese Vorgehensweise mit einem Effizienzverlust verbunden. Desweiteren bleibt der Einfluss der Partikelrotation auf das Feststoffverhalten unber¨ ucksichtigt. Insbesondere in wandnahen Bereichen sind durch zahlreiche Partikel-Wand-Kollisionen erh¨ohte Rotationsgeschwindigkeiten der Partikel zu beobachten [Som96]. Diese sind die Ursache einer am Partikel angreifenden Querkraft, welche die Partikelbewegung entscheidend beeinflussen kann. Auf diese Problematik wird in Kapitel 3 genauer eingegangen. Der Vorteil des Zwei-Kontinua-Verfahrens besteht darin, dass aufgrund einer makroskopischen Betrachtung beider Phasen selbst bei hohen Volumenanteilen der Feststoffphase eine effiziente Berechnung m¨oglich ist. Die f¨ ur die Gasphase aufgestellten Gleichungen k¨onnen in einer ¨ahnlichen Weise f¨ ur die Feststoffphase formuliert werden. F¨ ur diesen Gleichungstyp existieren leistungsstarke, numerische Methoden, woraus die besondere Effizienz dieses Verfahrens resultiert. Das im Rahmen dieser Arbeit neu entwickelte Verfahren zur Berechnung von Gas-Feststoff-Str¨omungen schließt das Zwei-Kontinua-Verfahren ein. Dabei wird der Vorteil der numerischen Effizienz ausgenutzt. Die wesentlichen Schwachstellen des Zwei-Kontinua-Verfahrens werden durch die Kopplung mit dem im n¨achsten Abschnitt beschriebenen Einzelpartikel-Verfahren vermieden. Einzelpartikel-Verfahren Die diskrete Natur der Feststoffphase in einer Gas-Feststoff-Str¨omung legt den Gedanken nahe, eine individuelle Betrachtung jedes einzelnen Partikels vorzunehmen. Das Einzelpartikel-Verfahren, in der Literatur auch unter dem Namen Euler/Lagrange-Verfahren bekannt, tr¨agt dieser mikroskopischen Betrachtungsweise Rechnung, indem die Bahnlinien der einzelnen Partikel unter dem Einfluss der auf sie wirkenden Kr¨afte bestimmt werden. Hierzu ist f¨ ur jedes Partikel die so genannte Basset-Boussinesq-Oseen Gleichung unter Einhaltung gegebener Randbedingungen auszuwerten. Die Beschreibung von Wechselwirkungen, wie der Partikel-Partikel- und Partikel-Wand-Kollisionen, erfolgt ebenfalls auf mikroskopischer Ebene durch entsprechende Modelle. Die mittleren Eigenschaften der Feststoffphase werden nach der Berechnung einer Vielzahl von Bahnlinien als Ergebnis einer statistischen Mittelung gewonnen. An dieser Stelle erkennt man bereits den Nachteil dieses Verfahrens, der sich aus der Notwendigkeit der Berechnung einer großen Anzahl von Partikelbahnen ergibt, die f¨ ur eine statistisch zuverl¨assige Mittelwertbildung unerl¨asslich sind. Diesem Nachteil steht eine Reihe von Vorteilen gegen¨ uber, die dieses Verfahren in den vergangenen Jahren 10

2.2 Mathematische Berechnungsverfahren zu einem wichtigen Werkzeug f¨ ur die Simulation von Gas-Feststoff-Str¨omungen gemacht haben. Hierzu geh¨ort insbesondere die M¨oglichkeit einer detaillierten Beschreibung der physikalischen Vorg¨ange, denen die Partikelbewegung unterliegt. So ist die Ermittlung der Partikelrotationsgeschwindigkeit, die f¨ ur die Berechnung der Partikelquerkr¨afte (Magnusund Saffman-Kraft) ben¨otigt wird, problemlos m¨oglich. Wechselwirkungen wie beispielsweise Partikel-Partikel-Kollisionen k¨onnen pr¨azise mit den Gesetzen der Stoßmechanik beschrieben werden. Ein weiterer Vorteil ist der breite Einsatzbereich des EinzelpartikelVerfahrens. So ist die Untersuchung eines polydispersen Feststoffsystems m¨oglich, indem eine Partikelgr¨oßenverteilung ber¨ ucksichtigt wird. Aufgrund der detaillierten Betrachtung der Vorg¨ange, die sich an einem einzelnen Partikel abspielen, k¨onnen zus¨atzliche Einfl¨ usse, wie beispielsweise eine Verbrennungsreaktion in die Modellierung mit einbezogen werden. Auch die Formulierung von Randbedingungen f¨ ur die Str¨omungsgr¨oßen der Partikelphase an der Bewandung stellt im Gegensatz zum Zwei-Kontinua-Verfahren keine Schwierigkeit dar [Som96]. Die oben genannte Einschr¨ankung des Einzelpartikel-Verfahrens, der hohe Zeitbedarf, offenbart sich insbesondere bei der Berechnung von dichten Gas-Feststoff-Str¨omungen. Hier kann allein die Anzahl der vorliegenden Partikel schnell Dimensionen erreichen, welche selbst die Kapazit¨at leistungsstarker Computer u ucksichtigung weiterer ¨bersteigt. Die Ber¨ Vorg¨ange wie z.B. der Partikel-Partikel-Kollisionen oder der R¨ uckwirkung der Partikel auf die Gasphase kann eine Simulation derartiger Gas-Feststoff-Str¨omungen aufgrund der sehr hohen Berechnungszeit unvertretbar machen [CST98]. Das im Rahmen dieser Arbeit neu entwickelte Verfahren tr¨agt zur L¨osung dieser Problematik bei, indem es die Berechnung der Feststoffphase mit dem Einzelpartikel-Verfahren auf einen begrenzten Bereich beschr¨ankt, innerhalb dessen die Anwendung des Zwei-Kontinua-Verfahrens aufgrund der im vorigen Abschnitt genannten Einschr¨ankungen nicht in Frage kommt. Dieser Bereich zeichnet sich durch ein hohes Maß an Unregelm¨aßigkeit in der Bewegung der Partikel aus und ist vorwiegend in der N¨ahe fester Hindernisse sowie an der Geometriebegrenzung vorzufinden. Hier k¨onnen die Vorteile des Einzelpartikel-Verfahrens, insbesondere die detaillierte Beschreibung des Feststoffverhaltens, ausgenutzt werden. Zudem er¨offnet die M¨oglichkeit der Beschr¨ankung der Simulation auf einen kleinen Bereich, in dem das Einzelpartikel-Verfahren zum Einsatz kommt, eine echte Perspektive zur Berechnung von dichten Gas-Feststoff-Str¨omungen innerhalb einer angemessenen Zeit.

11

2 Simulation von Gas-Feststoff-Str¨omungen

12

3

Formulierung der Bilanzgleichungen

Das im Rahmen dieser Arbeit entwickelte Verfahren zur Simulation von Gas-FeststoffStr¨omungen ist ein Hybridverfahren. Es kombiniert das Zwei-Kontinua-Verfahren mit dem Einzelpartikel-Verfahren. W¨ahrend die Beschreibung der Gasphase f¨ ur beide Verfahren identisch ist, existieren, wie in Kapitel 2 ausgef¨ uhrt, bei der Behandlung der Feststoffphase gravierende Unterschiede. In diesem Kapitel werden zun¨achst die zugrunde liegenden Bilanzgleichungen zur Berechnung der Gasphase vorgestellt. Anschließend werden ausgehend von der Formulierung der Boltzmann-Gleichung die beschreibenden Gleichungen der Feststoffphase nach dem Zwei-Kontinua- und dem Einzelpartikel-Verfahren hergeleitet. Dabei erfolgt die Formulierung der Bilanzgleichungen unter der Annahme isothermer Bedingungen und unter Ausschluss chemischer Reaktionen. Die Feststoffphase besteht aus ideal sph¨arischen, monodispersen Partikeln.

3.1 3.1.1

Beschreibung der Gasphase Bewegungsgleichungen

Die Bewegung der Gasphase wird durch die Bilanzgleichungen f¨ ur Masse, Impuls und Energie beschrieben. Dar¨ uber hinaus ist ein rheologisches Stoffgesetz sowie eine thermodynamische Zustandsgleichung zur vollst¨andigen Modellierung erforderlich. Betrachtet man ein beliebig geformtes, ortsfestes Kontrollvolumen Ω mit der Oberfl¨ache S, das von einer Gasphase mit der Geschwindigkeit ug durchstr¨omt wird (vgl. Abbildung 3.1), so ergibt sich f¨ ur die Massenbilanz der folgende Zusammenhang:   ∂ φg ρg dΩ + φg ρg ug · n dS = 0, (3.1) ∂t Ω S wobei dS ein infinitesimal kleines Element der Oberfl¨ache S und n den zu dS geh¨orenden nach außen gerichteten Normalenvektor darstellen. Der erste Term auf der linken Seite

S

n

Ω

ug

Abbildung 3.1: Allgemeines Kontrollvolumen. 13

3 Formulierung der Bilanzgleichungen von (3.1) beschreibt die Massen¨anderung im Kontrollvolumen Ω innerhalb der Zeitspanne ∂t, der zweite Term kennzeichnet den konvektiven Transport von Masse durch die Begrenzungsfl¨ache S. Da in den hier untersuchten Str¨omungskonfigurationen keine Phasen¨ uberg¨ange ber¨ ucksichtigt werden, finden sich auf der rechten Seite von (3.1) keine Quelloder Senkenterme. Geht man nun von einer Inkompressibilit¨at der Gasphase aus, was bei niedrigen Str¨omungsgeschwindigkeiten (Ma < 0.2) physikalisch gerechtfertigt und numerisch erforderlich ist, l¨asst sich die Gasdichte ρg aus (3.1) eliminieren [Str91]. Wird nun f¨ ur das Kontrollvolumen Ω die Bilanzgleichung f¨ ur den Impuls formuliert, so erh¨alt man:      ∂ φg ρg ug dΩ + φg ρg ug ug · n dS = Tg · n dS + φg ρg g dΩ + kgp dΩ. (3.2) ∂t Ω S S Ω Ω ¨ Die linke Seite von (3.2) beschreibt die zeitliche Anderung des Impulses sowie den konvektiven Impulsfluss u ¨ber die Oberfl¨ache von Ω. Der erste Term der rechten Seite repr¨asentiert den diffusiven Impulsfluss u ¨ber die Oberfl¨ache. Der Spannungstensor der Gasphase Tg auf der rechten Seite von (3.2) gibt die Normal- und Tangentialspannungen an der Oberfl¨ache des Kontrollvolumens wieder und ist f¨ ur inkompressible, newtonsche Fluide durch das zugeh¨orige rheologische Stoffgesetz gegeben:       2 (3.3) Tg = − pg + µg − ζg ∇ · ug I + µg ∇ug + (∇ug )T . 3 Hierbei bezeichnet ζg die Volumenviskosit¨at, die f¨ ur kompressible Medien ber¨ ucksichtigt werden muss. F¨ ur das hier betrachtete inkompressible Fluid nimmt diese folglich den Wert Null an. Die Gr¨oßen pg und I kennzeichnen den statischen Druck innerhalb der Gasphase sowie den Einheitstensor. Das zweite und dritte Integral der rechten Seite von (3.2) beschreiben den Einfluss der Schwerkraft sowie einen Impulsaustauschterm, der sich aus der Wechselwirkung mit der dispersen Phase ergibt. Die Gr¨oßen g und kgp bezeichnen in diesem Zusammenhang die Erdbeschleunigung und eine Volumenkraft, die sich aus der Partikelwechselwirkung ergibt. Die Diskussion von kgp erfolgt in Abschnitt 3.2. Die oben formulierten Gleichungen f¨ ur Masse und Impuls sind f¨ ur eine Zweiphasenstr¨omung definiert. Betrachtet man eine einphasige Str¨omung (φg =1 und kgp=0), so reduzieren sich (3.1) und (3.2) zur Kontinuit¨atsgleichung und Navier-Stokes-Gleichung in der allgemein bekannten Form. Aufgrund isothermer Bedingungen ist eine L¨osung der Energiebilanz nicht erforderlich.

3.1.2

Modellierung der Turbulenz

Eine turbulente Str¨omung ist unter anderem dadurch charakterisiert, dass sich die Schwankungsbewegungen auf unterschiedlichen L¨angenskalen vollziehen. Die Abmessungen der einzelnen Turbulenzelemente k¨onnen in der Gr¨oßenordnung der Str¨omungsgeometrie liegen, bis hin zu den kleinsten Elementen, deren Ausdehnung durch die so genannte Kolmogorovsche L¨angenskala η gegeben ist. Obwohl die kleinen Turbulenzelemente zahlenm¨aßig 14

3.1 Beschreibung der Gasphase u ¨berwiegen, ist der gr¨oßte Teil der Turbulenzenergie in den großskaligen Wirbeln enthalten [Fr¨o00]. Diese werden von der mittleren Str¨omung erzeugt und zerfallen in immer kleinere Wirbel. Mit abnehmender Wirbelgr¨oße werden viskose Kr¨afte dominant, welche die kinetische Energie der Turbulenzbewegung irreversibel in innere Energie der Str¨omung u uhren. Der Zusammenhang zwischen dem Energieinhalt und der Wirbelgr¨oße wird ¨berf¨ auch als Kolmogorov’sche Energiekaskade bezeichnet [Leo74]. Die in Abschnitt 3.1.1 formulierten Transportgleichungen gelten f¨ ur laminare und turbulente Str¨omungsvorg¨ange. Eine vollst¨andige ¨ortliche Beschreibung des turbulenten Str¨omungsfeldes setzt ein Rechengitter voraus, das auch noch die kleinsten Wirbel hinreichend genau erfasst. Dar¨ uber hinaus sollte der Zeitschritt auf eine Gr¨oße reduziert werden, welche das turbulente Zeitmaß aufl¨osen kann. Ein Verfahren, welches die erforderliche ¨ortliche und zeitliche Aufl¨osung besitzt, um alle Turbulenzeffekte erfassen zu k¨onnen, wird als die Direkte Numerische Simulation (DNS) bezeichnet. Hierbei werden die dreidimensionalen, zeitabh¨angigen Transportgleichungen durch eine entsprechende Diskretisierung des Str¨omungsfeldes in den Kontrollvolumina direkt d.h. ohne den Einsatz eines Turbulenzmodells gel¨ost. Die Ausdehnung der Kontrollvolumina h¨angt von der Gr¨oße der kleinsten turbulenten Wirbel ab, die durch die Kolmogorovsche L¨angenskala gegeben ist. Bei gr¨oßer werdender Turbulenzintensit¨at der Str¨omung nimmt der Wert dieses L¨angenmaßes ab. Dadurch erh¨oht sich der Rechenaufwand drastisch, da eine h¨ohere Aufl¨osung des Str¨omungsgebietes erforderlich ist. Aus diesem Grund ist die Anwendung der Direkten Numerischen Simulation bislang auf einfache Str¨omungsgeometrien (z.B. Kanalstr¨omung) bei moderater Turbulenzintensit¨at beschr¨ankt und somit f¨ ur die Praxis zur Zeit noch unbedeutend. Ein Verfahren, das den Rechenaufwand im Vergleich zu einer DNS erheblich reduziert, ist die Large-Eddy-Simulation (LES). Hierbei werden nur die großen, energiereichen Wirbel aufgel¨ost, w¨ahrend der Einfluss kleiner turbulenter Strukturen durch entsprechende Modelle (Feinstrukturmodelle) ber¨ ucksichtigt wird [Sma63, Fer96, DS97]. Ein bekanntes Feinstrukturmodell ist das Smagorinsky-Modell, das f¨ ur den Fall einer isotropen Turbulenz ¨ eine gute Ubereinstimmung mit experimentellen Untersuchungen liefert [BFR80]. Trotz eines bereits erheblich geringeren Rechenaufwandes der LES gegen¨ uber einer DNS, ist die Large-Eddy-Simulation derzeit noch als ein Verfahren zur Berechnung vergleichsweise einfacher Str¨omungskonfigurationen zu sehen. Ihr großes Potential bei der Analyse zahlreicher Str¨omungsph¨anomene l¨asst jedoch vermuten, dass sich die LES in naher Zukunft zu einem wichtigen Werkzeug f¨ ur die Simulation praxisrelevanter, turbulenter Str¨omungen entwickeln wird. In turbulenten Str¨omungen sind s¨amtliche Feldgr¨oßen Funktionen des Ortes und der Zeit. In vielen F¨allen, insbesondere bei industriellen Anwendungen, ist man lediglich an den Mittelwerten dieser Gr¨oßen interessiert, da meist auch nur diese messtechnisch erfasst werden k¨onnen. Zur mathematischen Behandlung werden alle zeitabh¨angigen Gr¨oßen Φ in einen Mittelwert Φ und einen Schwankungswert Φ
aufgespalten: Φ(t, x) = Φ(x) + Φ
(t, x).

(3.4)

Der Mittelwert in (3.4) wird u uber ¨ber die Zeit T bestimmt, die vergleichsweise ”lang” gegen¨ 15

3 Formulierung der Bilanzgleichungen der Zeit sein muss, in der sich die turbulenten Schwankungen vollziehen: Φ(x) =

1 T

T Φ(t, x) dt.

(3.5)

0

Wird nun f¨ ur die Feldgr¨oßen in (3.1) und (3.2) eine Aufspaltung gem¨aß (3.5) vorgenommen und eine zeitliche Mittelung der gesamten Gleichung durchgef¨ uhrt, so ergibt sich die so genannte Reynolds-gemittelte Massenbilanz:   ∂ φg dΩ + φg ug · n dS = 0, (3.6) ∂t Ω S wobei alle Korrelationen mit φ
g wegen ihrer geringen Bedeutung vernachl¨assigt werden. F¨ ur die Reynolds-gemittelte Impulsbilanz gilt entsprechend ∂ ∂t



 Ω

 φg ρg ug ug · n dS + φg ρg u
g u
g T · n dS S S   = Tg · n dS + φg ρg g dΩ + kgp dΩ,

φg ρg ug dΩ +

S

Ω

mit Tg als Reynolds-gemitteltem Spannungstensor:     2 Tg = − pg + µg ∇ · ug I + µg ∇ug + (∇ug )T . 3

(3.7)

Ω

(3.8)

Die Reynolds-gemittelte Impulsbilanz (3.7) unterscheidet sich von der Impulsbilanz (3.2) nur durch die Korrelation der Schwankungsgr¨oßen TR = −ρg u
g u
g T , die als Reynoldsspannungen bezeichnet werden. Sie sind durch die Mittelung des nichtlinearen Konvektionsterms entstanden und k¨onnen als zus¨atzliche Scherspannungen interpretiert werden, die durch die turbulenten Geschwindigkeitsfluktuationen hervorgerufen werden. Die Berechnung dieser Korrelationen erfordert weitere Informationen, die aus Turbulenzmodellen gewonnen werden k¨onnen. Dabei kann die Formulierung geeigneter Modellans¨atze auf zwei unterschiedlichen Wegen erfolgen. Zum einen lassen sich so genannte ReynoldsspannungsModelle herleiten, die f¨ ur die Bestimmung der Schwankungsgr¨oßen exakte Gleichungen bereitstellen [LRR75, Rodi87]. Infolgedessen zeichnet sich dieser Ansatz durch eine hohe Genauigkeit aus, ist jedoch aufgrund seiner Komplexit¨at (mehr als 6 zus¨atzliche partielle Differential-Gleichungen mit Termen h¨oherer Ordnung) rechenaufw¨andig und anspruchsvoll bez¨ uglich der Definition von Randbedingungen. Zum anderen k¨onnen die unbekannten Korrelationen u ¨ber einen Gradientenansatz ermittelt werden. Der Vorteil dieser Methode liegt darin, dass die Koeffizienten empirisch ermittelt oder aus halbempirischen Bestimmungsgleichungen berechnet werden k¨onnen, wobei die Anzahl zus¨atzlicher Transportgleichungen im Allgemeinen gering ist. Grundgedanke des Gradientenansatzes ist die Wirbelviskosit¨atshypothese nach Bousinesque, wonach die Reynoldsspannungen TR aus den 16

3.1 Beschreibung der Gasphase zeitlich gemittelten Geschwindigkeitsgradienten abgeleitet werden k¨onnen [LS72]:   2 TR = −ρg u
g u
g T = µt ∇ug + (∇ug )T − ρg k I. 3

(3.9)

Die turbulente kinetische Energie k ist hierbei definiert als: 1 k = (u
g u
g + vg
vg
+ wg
wg
). 2

(3.10)

Die Gr¨oße µt kennzeichnet die so genannte turbulente Viskosit¨at, die im Gegensatz zu der molekularen Viskosit¨at µg keine Materialkonstante, sondern eine ¨ortlich und zeitlich variable Gr¨oße darstellt. F¨ ur deren Berechnung existieren unterschiedliche Ans¨atze. Je nach der Anzahl der zus¨atzlichen Transportgleichungen werden diese in Null-, Ein- und Zweigleichungsmodelle unterteilt. Im Rahmen dieser Arbeit wird das von Launder & Spalding entwickelte und in der industriellen Praxis weit verbreitete k-ε-Modell, ein Zweigleichungsmodell, verwendet [LS72]. Hierzu werden die Transportgleichungen f¨ ur die turbulente kinetische Energie k und die Dissipationsrate ε gel¨ost. Die Dissipationsrate ε ist dabei wie folgt definiert: ε=

µg ∇u
g : (∇u
g )T . ρg

(3.11)

Mit diesen Gr¨oßen l¨asst sich die turbulente Viskosit¨at nach einer Kolmogorov-PrandtlBeziehung ermitteln [Kol42, Pra45]: µT = Cµ ρg

k2 , ε

(3.12)

ur die turbulente wobei Cµ eine empirische Konstante darstellt. Die Transportgleichungen f¨ kinetische Energie und die Dissipationsrate werden aus den gemittelten und zeitabh¨angigen Impulserhaltungss¨atzen hergeleitet [Hin75]. Unbekannte Korrelationen m¨ ussen dann noch durch die bekannten Gr¨oßen ug , k und ε ersetzt werden. Als Ergebnis erh¨alt man die Transportgleichungen f¨ ur die turbulente kinetische Energie k     ∂ µt ρg k dΩ + ρg k ug · n dS = ∇k · n dS + TR : ∇ug − ρg ε dΩ (3.13) ∂t Ω S S σk Ω und f¨ ur die Dissipationsrate ε     ∂ µt ε ε2 ρg ε dΩ+ ρg ε ug ·n dS = ∇ε·n dS + Cε1 TR : ∇ug −Cε2 ρg dΩ, (3.14) ∂t Ω σ k k S S ε Ω wobei der Ausdruck TR : ∇ug das Skalarprodukt zweier Tensoren kennzeichnet. Die Konstanten in (3.12), (3.13) und (3.14) wurden durch Regression an eine Vielzahl experimenteller Ergebnisse angepasst. Ein h¨aufig verwendeter Konstantensatz lautet [LS72]: 17

3 Formulierung der Bilanzgleichungen Cµ 0.09

Cε1 1.44

Cε2 1.92

σk 1.0

σε 1.3

Das Standard-k-ε Modell (3.13, 3.14) reduziert die Berechnung der komplexen Vorg¨ange der Turbulenz auf die L¨osung zweier Bilanzgleichungen. Aufgrund der getroffenen Annahmen (z.B. Isotropie der Turbulenz) und des Einsatzes von empirisch angepassten Konstanten sind dem Modell Grenzen hinsichtlich seiner G¨ ultigkeit gesetzt. Eine Anpassung ¨ ist zum Beispiel bei der Anwesenheit einer zweiten Phase erforderlich. Uber den Einfluss der dispersen Phase auf die Turbulenz gibt es jedoch widerspr¨ uchliche Aussagen, die eine Modellerweiterung erschweren. Es existieren zwar zahlreiche Vorschl¨age das k-ε-Modell zu erweitern, derartige Korrekturen sind aber nur von eingeschr¨anktem Nutzen, da sie f¨ ur spezielle Anwendungsf¨alle entwickelt wurden und somit keine Allgemeing¨ ultigkeit besitzen [CW86, TF94, GNN83]. Trotz dieser Einschr¨ankungen stellt das k-ε-Modell einen guten Kompromiss zwischen der erreichbaren Genauigkeit und dem erforderlichen Berechnungsaufwand dar. Dar¨ uber hinaus ergibt sich ein Vorteil bei der Anwendung im Zusammenhang mit mehrphasigen Str¨omungen, da die Modellkonstanten bereits gut angepasst sind und f¨ ur verschiedene Anwendungsf¨alle validiert wurden.

3.2

Beschreibung der Feststoffphase

Die Beschreibung des Bewegungsverhaltens der dispersen Feststoffphase erfolgt in Analogie zur Beschreibung der Molek¨ ulbewegung in einem idealen Gas. Hier wird die Annahme getroffen, die Molek¨ ule stellen harte Kugeln dar mit ausschließlich translatorischen Freiheitsgraden der Bewegung. Die intermolekularen Kollisionen sind bin¨ar und vollst¨andig elastisch. Der Bewegungszustand eines Molek¨ ulkollektivs l¨asst sich allgemein durch die Angabe von Ort und Geschwindigkeit jedes einzelnen Teilchens zu jedem Zeitpunkt beschreiben. Eine Funktion, die den Zustand eines Gases vollst¨andig beschreibt, ist die Wahrscheinlichkeitsdichteverteilung f (x, v, t). Das Produkt f (x, v, t) dx dv

(3.15)

gibt die Anzahl der Molek¨ ule an, deren Schwerpunkt sich zum Zeitpunkt t im Volumen (x, x + dx) befindet und die eine Geschwindigkeit im Intervall (v, v + dv) aufweisen. Die zeitliche und r¨aumliche Evolution von f wird durch die Boltzmann-Gleichung beschrieben (z.B. nach [Gid94]):   ∂f ∂f ∂f ∂f +v· +F· = . (3.16) ∂t ∂x ∂v ∂t coll Dabei stellt die linke Seite von (3.16) die kinetischen Terme dar, welche den Transport von f aufgrund einer ungest¨orten Molek¨ ulbewegung sowie einer durch die ¨außeren Kr¨afte beeinflussten Bewegung beschreiben. Die Summe der ¨außeren Kr¨afte wird in F zusammengefasst. Die rechte Seite von (3.16) repr¨asentiert die durch bin¨are St¨oße zwischen den ¨ Molek¨ ulen bedingte Anderung der Wahrscheinlichkeitsdichteverteilung f . 18

3.2 Beschreibung der Feststoffphase Eine analytische L¨osung der Boltzmann-Gleichung ist bislang nur unter vereinfachenden Annahmen m¨oglich. Ein Spezialfall stellt hierbei die Annahme eines thermodynamischen Gleichgewichts dar. Befindet sich ein Molek¨ ulkollektiv in einem thermodynamischen Gleichgewicht, d.h. es finden keinerlei Ausgleichsprozesse mehr statt, dann werden die drei Summanden auf der linken Seite von (3.16) zu Null. Damit (3.16) weiterhin erf¨ ullt ist, muss f¨ ur den als Kollisionsoperator bekannten Term auf der rechten Seite dieser Gleichung gelten:   ∂f ≡ 0. (3.17) ∂t coll Dies ist dann der Fall, wenn die Wahrscheinlichkeitsdichtefunktion f die Form der bekannten Maxwell-Verteilung fM einnimmt:   ρ |v − u|2 fM (v) = . (3.18) − 3 exp 2RT (2πRT ) 2 Dabei beschreiben ρ die Gasdichte, u die mittlere Molek¨ ulgeschwindigkeit und T die thermodynamische Temperatur, die ein Maß f¨ ur die kinetische Energie der Fluktuationsbewegung darstellt. Die Maxwell-Verteilung stellt somit die L¨osung der Boltzmann-Gleichung f¨ ur den Fall eines thermodynamischen Gleichgewichts dar. Unter dieser Voraussetzung lassen sich nun aus (3.16) Bilanzgleichungen herleiten, die den Transport von Masse, Impuls und Energie in einem idealen Gas beschreiben. Die hier kurz skizzierte Methodik der Beschreibung des Bewegungsverhaltens idealer Gase geh¨ort zum Gebiet der kinetischen Gastheorie. Eine ausf¨ uhrliche Beschreibung dieser Theorie ist dem Standardwerk von Chapman & Cowling zu entnehmen [CC70]. Um den Bewegungszustand eines dispersen Feststoffs mit vergleichbaren Methoden zu beschreiben, wurde die kinetische Gastheorie zu Beginn der achtziger Jahre zur kinetischen Theorie granularer Fluide (KTGF) erweitert. In dem folgenden Abschnitt 3.2.1 werden die Grundz¨ uge dieser Methode sowie die sich hieraus ergebenden Transportgleichungen f¨ ur die Feststoffphase vorgestellt. Eine andere M¨oglichkeit der L¨osung von (3.16) besteht in der Verwendung stochastischer Verfahren. Hierbei wird die zeitliche Evolution eines Molek¨ ulkollektivs mit stochastischen Methoden approximiert. Die Beschreibung der dazu notwendigen Techniken sowie der zu ber¨ ucksichtigenden Annahmen, insbesondere im Hinblick auf disperse Partikelkollektive, sind Gegenstand von Abschnitt 3.2.2.

3.2.1

Zwei-Kontinua-Verfahren

Beim Zwei-Kontinua-Verfahren wird angenommen, dass sowohl die Gas- als auch die Feststoffphase als zwei koexistierende Kontinua dargestellt werden k¨onnen. Die zeitlich gemittelten Erhaltungsgleichungen der Gasphase (3.6, 3.7) sind bereits in Abschnitt 3.1.2 vorgestellt worden. Die Transportgleichungen der Feststoffphase ergeben sich aus einer zeitlichen Mittelung der so genannten Maxwellschen Transportgleichungen, die ihren Ursprung in 19

3 Formulierung der Bilanzgleichungen der Boltzmann-Gleichung haben. Die Herleitung dieser Beziehungen erfolgt im Rahmen der kinetischen Theorie granularer Fluide, die im Gegensatz zur kinetischen Gastheorie eine Betrachtung von dispersen Feststoffsystemen erlaubt. Da die KTGF ausf¨ uhrlich von Gidaspow beschrieben wurde [Gid94], sollen an dieser Stelle lediglich die zu treffenden Annahmen sowie eine Skizze der Herleitung bis hin zur Formulierung der Transportgleichungen der Feststoffphase gegeben werden. Eine zweiphasige Gas-Feststoff-Str¨omung unterscheidet sich von der Str¨omung eines ¨ Molek¨ ulkollektivs durch die Anwesenheit einer zweiten Phase. Daher sind trotz der Ahnlichkeit des Bewegungsverhaltens des fluidisierten Feststoffs zu dem eines idealen Gases wesentliche Unterschiede festzustellen. Die wichtigsten lassen sich wie folgt zusammenfassen: • die Partikel-Partikel-Kollisionen erfolgen inelastisch, das heißt, um eine Fluidisierung der Feststoffphase aufrechtzuerhalten, muss der Verlust an kinetischer Energie infolge der Kollisionen durch einen Energietransfer von der Gasphase zur Feststoffphase kompensiert werden, • zwischen den Kollisionen wirken auf die Partikel ¨außere Kr¨afte, die zum einen aus der Druckverteilung um das Partikel, zum anderen aus den Schubspannungen an der Oberfl¨ache des Partikels resultieren, • die Gravitationskraft sowie die Rotationsbewegung der Partikel ist zu ber¨ ucksichtigen, • die turbulenten Strukturen der Gasphase beeinflussen die Fluktuationsbewegung der Partikel. Die hier aufgef¨ uhrten Unterschiede zwischen einem dispersen Feststoffsystem und einem Molek¨ ulkollektiv werden bis auf den Schwerkrafteinfluß und die Rotationsbewegung in der kinetischen Theorie granularer Fluide ber¨ ucksichtigt. Die Abweichungen vom thermodynamischen Gleichgewicht, die aus der Existenz treibender Kr¨afte herr¨ uhren, werden allgemein durch folgenden Ansatz erfasst [Gid94]: f = f (0) + f (1) = f (0) (1 + Φ(1) ).

(3.19)

f stellt die zweite Approximation zur Boltzmann-Gleichung dar (die Maxwell-Verteilung f (0) wird in diesem Zusammenhang auch als erste Approximation bezeichnet); der Ausdruck Φ(1) steht f¨ ur ein Korrekturpolynom, das bestimmte Abweichungen vom thermodynamischen Gleichgewicht ber¨ ucksichtigt. Die Herleitung sowie die Diskussion dieser Gr¨oße erfolgt in Kapitel 5. Unter Verwendung von (3.19) k¨onnen nun die Transportgleichungen der Feststoffphase in analoger Weise zum Vorgehen in der kinetischen Gastheorie hergeleitet werden. Hierzu wird die Boltzmann-Gleichung mit einer Gr¨oße q, die den Transport einer Str¨omungseigenschaft beschreibt, multipliziert und u ¨ber den gesamten Geschwindigkeitsraum integriert:       ∂f ∂f ∂f ∂f q +v· +F· dv = q dv. (3.20) ∂t ∂x ∂v ∂t coll 20

3.2 Beschreibung der Feststoffphase Das Ergebnis der Integration stellt die Maxwellschen Transportgleichungen dar. F¨ ur q = φp ρp ergibt sich die Massenerhaltungsgleichung der Partikelphase:   ∂ φp dΩ + φp up · n dS = 0, (3.21) ∂t Ω S mit dem mittleren Volumenanteil φp und der mittleren Geschwindigkeit der Partikel up . Setzt man nun q = φp ρp up , so kann die Impulserhaltungsgleichung wie folgt angegeben werden:      ∂ φp ρp up dΩ + φp ρp up up · n dS = Tp · n dS + φp ρp g dΩ − kgp dΩ, (3.22) ∂t Ω S S Ω Ω wobei Tp den Spannungstensor der Partikelphase darstellt:       2 µp − ζp ∇ · up I + µp ∇up + (∇up )T . Tp = − pp + 3

(3.23)

Hier bezeichnen pp , ζp und µp den Partikeldruck, die Volumen- und die Scherviskosit¨at der Partikelphase. Die drei Unbekannten sowie der Impulsaustauschterm kgp in (3.22), kgp = β (up − ug ), der als Produkt aus einer Widerstandsfunktion β und der Relativgeschwindigkeit zwischen Partikelensemble und Gas modelliert wird, sind weiter unten erl¨autert. Die dritte Maxwellsche Transportgleichung wird f¨ ur die Fluktuationsenergie der Partikelphase formuliert. Dabei wird angenommen, dass die Fluktuationsbewegungen der Partikel isotrop sind und sich deshalb aus der Definition der kinetischen Energie Ep ein Zusammenhang zwischen der Fluktuationsbewegung u
p und der granularer Temperatur Θ ergibt: 1 1 1 3 2 Ep = mp up 2 + mp u
p = mp up 2 + mp Θ. 2 2 2 2

(3.24)

Setzt man nun in (3.20) q = 3/2mp Θ, so kann folgende Bilanzgleichung f¨ ur die granulare Temperatur Θ gefunden werden:   ∂ 3 3 φp ρp Θ dΩ + φp ρp up · nΘ dS (3.25) ∂t Ω 2 S 2     = Tp : ∇up dΩ + κp ∇Θ · n dS − γp dΩ + Ψp dΩ. Ω

S

Ω

Ω

¨ Die linke Seite von (3.25) beschreibt die zeitliche Anderung und den konvektiven Transport der Fluktuationsenergie der Partikel. Die ersten beiden Terme der rechten Seite stellen die Produktion sowie den diffusiven Transport der granularen Temperatur dar. Die verbleibenden beiden Summanden bezeichnen die Dissipation sowie den Austausch der Fluktuationsenergie zwischen Gas und Partikel. Die Geschwindigkeitsgradienten der Partikelstr¨omung werden durch den Spannungstensor Tp (3.23) charakterisiert. Die hierin enthaltene Gr¨oße pp bezeichnet den Partikeldruck. 21

3 Formulierung der Bilanzgleichungen Er repr¨asentiert das Bestreben der Partikel sich innerhalb der Gasphase r¨aumlich auszubreiten und w¨ urde bei Abwesenheit anderer Kr¨afte zu einer gleichm¨aßigen Partikelbeladung im gesamten Str¨omungsgebiet f¨ uhren. Die Berechnung des Partikeldruckes erfolgt aus einem Kinetik- und einem Kollisionsterm [LSJ84]: pp = pp,kin + pp,koll .

(3.26)

Der kinetische Term erfasst die translatorische Bewegung der Partikel und ergibt sich zu: pp,kin = φp ρp Θ,

(3.27)

der Kollisionsterm ber¨ ucksichtigt Partikelkollisionen und lautet: pp,koll = 2(1 + e)g0φp 2 ρp Θ.

(3.28)

Hierbei bezeichnen e und g0 die Stoßzahl und die radiale Verteilungsfunktion. Die Stoßzahl e beschreibt den Verlust an kinetischer Energie bei einer Kollision. Ein vollkommen elastischer Stoß wird mit e = 1 beschrieben; e = 0 f¨ uhrt zu einer maximalen Dissipation kinetischer Energie w¨ahrend einer Kollision. Im Rahmen dieser Arbeit wird e als eine Materialkonstante behandelt, die Abh¨angigkeit von anderen Randbedingungen wie Auftreffwinkel oder Auftreffgeschwindigkeit wird vernachl¨assigt. Die radiale Verteilungsfunktion g0 ist ein Maß f¨ ur die Wahrscheinlichkeit, dass sich zwei Teilchen eines dichten Partikelkollektivs ber¨ uhren. F¨ ur die Modellierung dieser Gr¨oße existiert eine Vielzahl von ¨ Ans¨atzen; einen Uberblick dar¨ uber gibt Boemer [Boe96]. In dieser Arbeit wird der Ansatz von Ma & Ahmadi verwendet, der f¨ ur Suspensionen mit Volumenanteilen im Bereich von φp < 10−2 eine gute N¨aherung darstellt [MA90, Ind01]:  3 −0.678 

φp . g0 = 1 + 4φp 1 + 2.5φp + 4.59φp 2 + 4.515φp 3 1 − 0.644

(3.29)

Eine weitere Unbekannte in (3.23) ist die Scherviskosit¨at µp , die keine Stoffgr¨oße ist, son¨ dern eine Variable, die vom Str¨omungszustand abh¨angig ist. Ahnlich wie beim Partikeldruck setzt sie sich aus einem kinetischen und einem kollisionalen Term zusammen. Das Ergebnis kann wie folgt angegeben werden [Gid94]:
  2 Θ 10 π 4 4 2 µp = ρp dp (3.30) 1 + (1 + e)g0 φp + (1 + e)g0φp . π 96 (1 + e)g0 5 5 Die Volumenviskosit¨at ζp resultiert aus einer lokalen Kompression der Partikelphase und wird in folgender Form verwendet [LSJ84]:
Θ 4 ζp = φp 2 ρp dp g0 (1 + e) . (3.31) 3 π 22

3.2 Beschreibung der Feststoffphase Der zweite Term der rechten Seite von (3.25) stellt den diffusiven Transport der granularen Temperatur dar. Die darin enthaltene Gr¨oße κp wird als W¨armeleitf¨ahigkeit der Partikelphase bezeichnet und kann angegeben werden als [Gid94]:
  2 Θ 150 π 6 (3.32) κp = ρp dp 1 + (1 + e)g0 φp + 2φp 2 g0 (1 + e) . π 384 (1 + e)g0 5 Die Dissipationsrate γp in (3.25) beschreibt die Energiedissipation aufgrund inelastischer Partikelkollisionen [Gid94]: 
4 Θ − ∇ · up . (3.33) γp = 3(1 − e2 )φp 2 ρp g0 Θ dp π Bei v¨ollig elastischen St¨oßen (e = 1) findet keine Dissipation kinetischer Energie statt ucksichtigt die Produktion (γp = 0). Die Gr¨oße Ψ im letzten Summanden in (3.25) ber¨ bzw. Dissipation der Fluktuationsenergie Θ aufgrund von Reibungskr¨aften zwischen den Partikeln und der Gasphase: Ψp = β(Kgp − 3Θ).

(3.34)

upfung Die Gr¨oße Kgp wird als Geschwindigkeitskovarianz bezeichnet und stellt eine Verkn¨ zwischen der Gas- und Partikelfluktuation dar. Sie wird in Kapitel 6 genauer beschrieben. Die im Rahmen des Zwei-Kontinua-Verfahrens noch nicht n¨aher spezifizierte Gr¨oße ist der Impulsaustauschterm kgp aus (3.7) und (3.22). Dieser beschreibt die Impuls¨anderung aufgrund der Wechselwirkung zwischen der Gas- und Partikelphase. Die in diesem enthaltenen wesentlichen Kraftbeitr¨age sind: • die Widerstandskraft FW und • die Querkr¨afte (Magnus- und Saffman-Kraft). Weitere Kr¨afte, wie die virtuelle Massenkraft oder die Basset-Kraft ber¨ ucksichtigen die instation¨are Bewegung der Partikel und den daraus resultierenden Einfluss des Mediums, welches das Partikel unmittelbar umgibt. Diese Kraftbeitr¨age k¨onnen jedoch aufgrund ¨ des hohen Dichteverh¨altnisses ρp /ρg vernachl¨assigt werden [KES98]. Ahnliches gilt f¨ ur die Druckkraft, die sich aus dem dynamischen Druckgradienten in der Gasstr¨omung ergibt. Auch diese kann aufgrund ihrer geringen Bedeutung in turbulenten Gas-Feststoff-Str¨omungen vernachl¨assigt werden [Som96]. Im Gegensatz dazu k¨onnen die Querkr¨afte einen entscheidenden Einfluss auf die Partikelbewegung haben. Deren Bedeutung h¨angt jeweils von den Randbedingungen der Str¨omung ab. Wichtig in diesem Zusammenhang ist vor allem die Geometrie der Str¨omung sowie die Feststoffeigenschaften wie z.B. der Partikeldurchmesser. Beide Querkr¨afte entstehen durch eine unterschiedliche Anstr¨omgeschwindigkeit an der Oberfl¨ache des Partikels, so dass sich eine unsymmetrische Druckverteilung auf der Partikeloberfl¨ache ausbildet. Die resultierende Kraft wirkt dann in Richtung des geringeren Druckes. Die Saffman-Kraft, die aus einem Geschwindigkeitsgradienten der Str¨omung 23

3 Formulierung der Bilanzgleichungen resultiert, kann insbesondere in Wandn¨ahe große Werte annehmen [Som96]. Die MagnusKraft ergibt sich aus der Rotation der Partikel im Str¨omungsfeld der Gasphase. Hohe Rotationsgeschwindigkeiten k¨onnen vor allem in solchen Str¨omungsgeometrien auftreten, wo begrenzende W¨ande oder Hindernisse h¨aufige Partikel-Wand-Kollisionen induzieren. Da im Rahmen dieser Arbeit die Umstr¨omung von Hindernissen in einer Kanalgeometrie untersucht wird, darf der Einfluss der Magnus- sowie der Saffman-Kraft in der N¨ahe dieser Begrenzungen auf keinen Fall vernachl¨assigt werden. Im Kernbereich der Str¨omung hingegen, wo die Partikelbewegung vorwiegend von aerodynamischen Kr¨aften dominiert wird und der Schergradient der Str¨omung gering ist, kann der Beitrag beider Querkr¨afte außer acht gelassen werden. Da das Zwei-Kontinua-Verfahren in dem zuletzt genannten Bereich der Str¨omung Anwendung findet (s. Kapitel 7), kann auf den Einsatz beider Kraftkomponenten innerhalb dieses Verfahrens verzichtet werden. Das Einzelpartikel-Verfahren, das in den u ¨brigen Bereichen der Str¨omungsgeometrie eingesetzt wird, muss dagegen den Einfluss beider Querkr¨afte auf das Partikelverhalten aus den genannten Gr¨ unden ber¨ ucksichtigen. Die Formulierung der beschreibenden Gleichungen dieser Kraftbeitr¨age erfolgt im folgenden Abschnitt zusammen mit der Beschreibung der Widerstandskraft, die in beiden Verfahren Ber¨ ucksichtigung findet.

3.2.2

Einzelpartikel-Verfahren

Wie bereits in Abschnitt 3.2 erw¨ahnt, ist man nur in Ausnahmef¨allen in der Lage eine analytische L¨osung der Boltzmann-Gleichung anzugeben. Um dennoch praxisrelevante Probleme, wie beispielsweise den Wiedereintritt eines Raumgleiters in die Erdatmosph¨are, behandeln zu k¨onnen, wurden stochastische Simulationsverfahren entwickelt. Diese wurden zun¨achst f¨ ur die Untersuchung von Molek¨ ulkollektiven eingesetzt und approximieren die zeitliche Evolution eines Partikelsystems [Bird76]. Das Einzelpartikel-Verfahren, das im Folgenden vorgestellt wird, erm¨oglicht die Beschreibung des Bewegungsverhaltens eines Partikelkollektivs in einer Gas-Feststoff-Str¨omung. Dabei wird auf Methoden sowie Simulationstechniken zur¨ uckgegriffen, die denen zur Beschreibung von Molek¨ ulsystemen weitestgehend entsprechen. Die Simulation der Bewegung eines Partikelkollektivs basiert auf einer von Bird vorgeschlagenen Splitting-Technik [Bird76]. In einem ersten Schritt erfolgt die simultane Berechnung einer Vielzahl von Partikelbahnen ohne Ber¨ ucksichtigung interpartikul¨arer Kollisionen: t+∆t 

xj (t + ∆t) = xj (t) +

vj (s)ds.

(3.35)

t

Im Zeitintervall ∆t, in dem die Berechnung stattfindet, ¨andert sich die Partikelgeschwindigkeit vj nur aufgrund der auf das Partikel wirkenden Kr¨afte oder der Wechselwirkung mit der Geometrieberandung. Dabei muss die Anzahl der simulierten Partikel groß genug sein, um statistisch zuverl¨assige Aussagen u ¨ber das mittlere Bewegungsverhalten der Partikelphase treffen zu k¨onnen. Der zweite Schritt beinhaltet die Berechnung der interparti24

3.2 Beschreibung der Feststoffphase kul¨aren Kollisionen, welche die Partikel im Verlaufe des vergangenen Zeitschrittes erfahren haben. Die Bestimmung des mittleren Bewegungsverhaltens der Partikel erfolgt durch eine Auswertung der, durch das Einzelpartikel-Verfahren bestimmten, Approximation der Verteilungsfunktion f (t, x, v). Translatorische Partikelbewegung ¨ Die Anderung des Bewegungszustandes eines Einzelpartikels resultiert aus den an diesem Partikel angreifenden Kr¨aften und Momenten. Betrachtet man die translatorische Bewegung eines Partikels in einer Gasstr¨omung, welches zun¨achst keinerlei Wechselwirkung mit ¨ anderen Partikeln oder der Berandung erf¨ahrt, so kann die zeitliche Anderung seines Impulses durch die folgende Gleichung beschrieben werden: mp

dvp = FG + FA + FW + FS + FM . dt

(3.36)

Die ersten beiden Kraftbeitr¨age auf der rechten Seite von (3.36) stellen die K¨orperkr¨afte Gewichtskraft FG und Auftriebskraft FA - dar: FG =

π ρp dp 3 g, 6

FA =

π ρg dp 3 g. 6

(3.37)

Die verbleibenden drei Kraftkomponenten werden als aerodynamische Kr¨afte bezeichnet; es sind die Widerstandskraft, die Saffman-Kraft und die Magnus-Kraft. Die Widerstandskraft FW , die auf ein Partikel in einer parallelen Fluidstr¨omung wirkt, resultiert aus einem Reibungs- und einem Druckwiderstand. Die Gleichung zur Bestimmung des Gesamtwiderstandes lautet: FW = cd

ρg πdp 2 |vrel |vrel , 2 4

vrel = ug − vp ,

(3.38)

wobei vrel die Relativgeschwindigkeit zwischen Fluid und Partikel ist. Der Widerstandsbeiwert cd ist eine Funktion der mit vrel gebildeten Partikel-Reynolds-Zahl Rep Rep =

ρp dp |vrel | µg

(3.39)

und kann mit Hilfe der Korrelation von Turton & Levenspiel f¨ ur Rep < 2·105 angegeben werden [TL86]: cd =

0.413 24 1 + 0.173 Rep0.657 + . Rep 1 + 16300 Re−1.09 p

(3.40)

Eine Modifikation dieser Gleichung ergibt sich aus der Wechselwirkung eines Partikels mit den turbulenten Strukturen der Gasphase. Es wird angenommen, dass der Einfluss der Fluidturbulenz auf den Widerstandsbeiwert durch Interaktionen der kleinen Wirbel mit 25

3 Formulierung der Bilanzgleichungen der Partikeloberfl¨ache zustande kommt [BGM98]. Dies f¨ uhrt zu einer Erh¨ohung des Widerstandsbeiwertes. Eine Beziehung die diesen Zusammenhang beschreibt lautet [BGM98]:   3 dp −4 cdt = cd 1 + 8.76 · 10 , (3.41) η mit η als der Kolmogorow´schen-L¨angenskala:  η=

µg 3 ρg 3 ε

 14 .

(3.42)

Ermittelt man nun in einer Volumeneinheit die mittlere Widerstandskraft FW sowie die mittlere Relativgeschwindigkeit vrel = ug − up , so kann die im letzten Abschnitt 3.2.1 eingef¨ uhrte Widerstandsfunktion β zur Ber¨ ucksichtigung der Impulswechselwirkung zwischen Gas und Partikel aus folgender Gleichung berechnet werden: β=

np |FW | , |vrel |

(3.43)

wobei np die lokale Anzahldichte darstellt. Die Bedeutung der Querkr¨afte, der Saffmanund der Magnus-Kraft, wurde bereits in Abschnitt 3.2.1 erl¨autert. Beide Kraftkomponenten werden in der Simulation mit dem Einzelpartikel-Verfahren ber¨ ucksichtigt. Die Gleichung zur Bestimmung der Saffman-Kraft lautet: FS = 6.46

d2p 4




ρg µg (vrel × ωg ) f (Rep , ReS ) , |ωg |

(3.44)

wobei ωg = ∇ × ug /2 der Vektor der Rotationsgeschwindigkeit der Gasphase ist. Der Korrekturfaktor f (Rep , ReS ) erlaubt die Anwendung dieser Gleichung auch bei h¨oheren Partikel-Reynolds-Zahlen. Nach Mei gilt [Mei92]: ⎧   √ √ Rep ⎨ ur Rep ≤ 40 (1 − 0.3314 ) exp − + 0.3314  f¨ f (Rep , ReS ) = , (3.45) 10  ⎩ f¨ ur Rep > 40 0.0524  Rep mit  = ReS /2Rep . ReS stellt die Reynolds-Zahl der Scherstr¨omung dar: ReS =

ρg d2p |ωrel | , µg

ωrel = ωg − ωp .

(3.46)

Die Gleichung zur Bestimmung der Magnus-Kraft lautet: FM = cM 26

π d2p ρg |vrel | (ωrel × vrel ) , 2 |ωrel | 4

(3.47)

3.2 Beschreibung der Feststoffphase wobei cM den Auftriebsbeiwert kennzeichnet. Die Bestimmung von cM ist in der Vergangenheit Gegenstand zahlreicher Untersuchungen gewesen. Die Ergebnisse sind unter anderem in der Arbeit von Sommerfeld verglichen worden und zeigen widerspr¨ uchliche Resultate [Som96]. In der vorliegenden Arbeit wird der Zusammenhang nach Tsuji et. al verwendet, da dieser den hier untersuchten Partikelreynoldszahlenbereich am besten beschreibt [TMM85]:  1 dp |ωrel | 0.4G falls G ≤ 1 cM = . (3.48) , G= 0.4 falls G > 1 2 |vrel | Der Einfluss weiterer Kr¨afte wie der virtuellen Massenkraft, der Basset-Kraft und der Druckkraft wird, wie schon beim Zwei-Kontinua-Verfahren (s. Abschnitt 3.2.1), vernachl¨assigt. Damit sind alle relevanten Kr¨afte erfasst, welche die translatorische Partikelbewegung entscheidend beeinflussen. Die resultierende Differentialgleichung f¨ ur die Partikelgeschwindigkeit lautet:  ∂vp 3 ρg |vrel | = (ωrel × vrel ) (3.49) cd |vrel |vrel + cM ∂t 4 ρp dp |ωrel |  0.5 2 µg + 6.46 cS (vrel × ωg ) + g. π ρg |ωg | Die L¨osung von (3.49) erfolgt mit Hilfe des Polygonzug-Verfahrens: ∆t  (t) F . (3.50) mp  (t) Hierbei bezeichnet F die Summe der im Zeitschritt t auf das Partikel wirkenden, a¨ußeren Kr¨afte, w¨ahrend ∆t den Zeitschritt der Berechnung repr¨asentiert. Die nach dem Zeitschritt ∆t neu eingenommene Partikelposition ergibt sich entsprechend zu: vp(t+∆t) ≈ vp(t) +

(t) ≈ x(t) x(t+∆t) p p + vp ∆t.

(3.51)

Rotatorische Partikelbewegung ¨ Die Anderung des Drehimpulses eines Partikels ist gleich der Summe aller an ihm angreifenden Momente Mi : dωp  = Mi . (3.52) Ip dt i Hierbei bezeichnet Ip das Tr¨agheitsmoment einer Kugel, f¨ ur das gilt: Ip =

1 mp dp 2 . 10

(3.53)

Die Summe der am Partikel angreifenden Momente resultiert aus der Bewegung des umgebenden Fluides sowie aus der Wechselwirkung mit W¨anden oder anderen Partikeln. 27

3 Formulierung der Bilanzgleichungen W¨ahrend der Einfluss der Kollisionen weiter unten beschrieben wird, l¨asst sich f¨ ur die Wirkung des Fluides auf die Rotation des Partikels eine von Dennis formulierte Beziehung angeben [DSI80]: M = cR

ρg dp 5 |ωrel |ωrel . 64

(3.54)

cR stellt den Rotationsbeiwert dar. Dessen Abh¨angigkeit von der Reynolds-Zahl der Rotation ReR = ρg dp 2

|ωrel | 4µg

l¨asst sich wie folgt zusammenfassen: ⎧ 16π ⎪ ⎪ f¨ ur 0 ≤ ReR < 101 ⎪ ⎪ ReR ⎪ ⎪ ⎪ ⎪ 6.45 32.1 ⎪ ⎪ ⎪ √ + f¨ ur 101 ≤ ReR < 103 ⎪ ⎪ Re Re R ⎪ R ⎪ ⎨ 6.8 √ f¨ ur 103 ≤ ReR < 4 · 104 . cR (ReR ) = ⎪ ReR ⎪ ⎪ ⎪ 0.058 ⎪ ⎪ ⎪ √ f¨ ur 4 · 104 ≤ ReR < 4 · 105 ⎪ 20 ⎪ ReR ⎪ ⎪ ⎪ ⎪ 0.397 ⎪ ⎪ f¨ ur 4 · 105 ≤ ReR < 107 ⎩ √ 5 ReR

(3.55)

(3.56)

Setzt man nun (3.54) in (3.52) ein, so ergibt sich folgende Differentialgleichung zur Bestimmung der Rotationsgeschwindigkeit eines Partikels: dωp 15 ρg cR = |ωrel |ωrel . dt 16 ρp π

(3.57)

Die L¨osung von (3.57) erfolgt wie im translatorischen Fall mit Hilfe des PolygonzugVerfahrens. Einfluss der Fluidturbulenz auf die Partikelbewegung Die Berechnung der Partikelbewegung in einer turbulenten Gas-Feststoff-Str¨omung erfordert die Kenntnis der Momentangeschwindigkeit der Gasphase am Aufenthaltsort des Partikels. Wie bereits in Abschnitt 3.1.2 erw¨ahnt, setzt sich diese aus einem zeitlichen Mittelwert und einem Schwankungsanteil zusammen. Die zeitlich gemittelte Gasgeschwindigkeit ergibt sich direkt aus der L¨osung der in Abschnitt 3.1 vorgestellten Gleichungen zur Berechnung der Gasphase; der Schwankungsanteil, der insbesondere bei kleinen Partikeln zu verst¨arkter Partikeldispersion f¨ uhren kann, muss durch ein geeignetes Modell erfasst werden. Im Rahmen dieser Arbeit wird ein stochastisches Modell, das so genannte Diskrete-Wirbel-Modell 28

3.2 Beschreibung der Feststoffphase verwendet [Mil90]. Hierin werden die Turbulenzeigenschaften der Gasphase entlang einer jeden Partikeltrajektorie ber¨ ucksichtigt. Die Verteilung der Schwankungsgeschwindigkeit an einem festen Ort wird unter der Annahme isotroper Turbulenz durch eine Gaussverteilung mit der Standardabweichung σg modelliert:
2 σg = k, (3.58) 3 wobei k die turbulente kinetische Energie der Gasphase ist. Es wird nun angenommen, dass die aus einem Zufallsprozeß erzeugte Schwankungsgeschwindigkeit, der ein Partikel ausgesetzt ist, u ¨ber einen festgelegten Zeitraum ∆tint entlang der Partikelbahn konstant bleibt. Dabei ist ∆tint durch das Minimum aus der Wirbelzerfallszeit TE und der Durchquerungszeit TD gegeben: ∆tint = min (TE , TD ) ,

(3.59)

mit TE = 0.3k/ε und TD = TE σg /|vrel |. Ist ∆tint gr¨oßer als eine der beiden Zeitskalen, so wird eine neue Schwankungsgeschwindigkeit bestimmt, mit der das Partikel in Wechselwirkung tritt. Das Diskrete-Wirbel-Modell stellt eine einfache und effiziente Methode dar, Fluidfluktuation entlang der Partikelbahn zu berechnen und ist in der Lage die Partikeldispersion selbst in komplexeren zweiphasigen Str¨omungen zu beschreiben. Einfluss der Partikel auf die Fluidbewegung ¨ Die Anderung der Geschwindigkeit eines Partikels aufgrund des Impulsaustausches mit dem umgebenden Fluid impliziert, dass auch die Bewegung der Gasphase durch die Bewegung der Partikel beeinflusst wird. Um dem Impuls¨ ubertrag bei der Berechnung der Gasphase Rechnung zu tragen, wurde von Crowe die so genannte Particle-Source-In-Cell (PSIC) Methode entwickelt [CSS77]. Hierin wird die Impuls¨anderung aller Partikel aufsummiert, die ein gegebenes Kontrollvolumen ∆Ω im Zeitschritt ∆t durchqueren. Als Ergebnis erh¨alt man den in (3.7) enthaltenen Impulsaustauschterm kgp: kgp =

Nτ N ∆Ω,j    1  (j) (j−1) − g ∆τj , mp,k vp,k − vp,k ∆Ω ∆t j=1

(3.60)

k=1

wobei ∆t die Summe der Subzeitschritte ∆τ repr¨asentiert, in denen die Berechnung der Partikelphase erfolgt. Die Impuls¨anderung in (3.60) wird um den Schwerkraftterm mp,k g∆τj verringert, da dieser zwar bei der Berechnung der Impuls¨anderung eines Partikels ber¨ ucksichtigt wird, jedoch kein Resultat der Gasphase ist. Partikel-Wand-Kollisionen Die Modellierung der Partikel-Wand-Kollisionen spielt bei der Beschreibung der Partikelbewegung in einer berandeten Gas-Feststoff-Str¨omung eine wichtige Rolle. So kann das 29

3 Formulierung der Bilanzgleichungen Bewegungsverhalten der Partikelphase durch vermehrte Wandkollisionen stark beeinflusst werden. Dies gilt insbesondere f¨ ur große Partikel, die aufgrund ihrer h¨oheren Tr¨agheit langsamer auf die Kr¨afte der Gasphase reagieren k¨onnen und somit h¨aufiger mit der Berandung kollidieren. Die dabei induzierte Partikelrotation hat durch die Wirkung der Magnus-Kraft einen entscheidenden Einfluss auf die Partikelbewegung. Der Kollisionsvorgang ist von mehreren Einflussparametern abh¨angig; die wichtigsten davon sind [Som96]: • der Auftreffwinkel, • die Translations- und Rotationsgeschwindigkeit des Partikels vor dem Aufprall, • die Materialpaarung zwischen Wand und Partikel, • die Form des Teilchens und • die Rauigkeit der Wand. Im Rahmen dieser Arbeit findet die Modellierung der Partikel-Wand-Kollisionen ohne Beschr¨ankung der Allgemeinheit an glatten W¨anden und unter Verwendung von sph¨arischen Partikeln statt. Die Berechnung der translatorischen und rotatorischen Partikelgeschwindigkeit nach einer Wandkollision erfolgt mit Hilfe der Impuls- und Drehimpulsgleichungen, die zu den bekannten Stoßgesetzen f¨ uhren. Hierzu werden Informationen u ¨ ber den Restitutionskoeffizienten eW und den Gleitreibungsbeiwert µW ben¨otigt. Beide Parameter sind im Wesentlichen von der Materialpaarung abh¨angig. Die in Kapitel 7 durchgef¨ uhrten Simulationen gehen von Polystyrol-Partikeln aus, die sich in einem Glaskanal bewegen. Aus Experimenten konnte hierf¨ ur ein Parametersatz von ew = 0.96 und µW = 0.201 gefunden werden [Ind01]. Eine ausf¨ uhrliche Beschreibung des Stoßvorganges kann der Arbeit von Sommerfeld entnommen werden [Som96]. Partikel-Partikel-Kollisionen Der Einfluss von Partikel-Partikel-Kollisionen auf das Bewegungsverhalten der Feststoffphase in einer dispersen Gas-Feststoff-Str¨omung offenbart sich bereits bei moderaten Beladungen. Experimentelle Untersuchungen von Tanaka et al. in einem vertikalen Rohr mit einem Feststoffvolumenanteil von φp = 4 · 10−4 haben gezeigt, dass interpartikul¨are Kollisionen zu einer verst¨arkten Dispersion der Partikel und damit zu einer gleichm¨aßigeren Konzentrationsverteilung im Rohrquerschnitt f¨ uhren [TTT89]. Diese Messungen wurden durch entsprechende Simulationen von Tanaka & Tsuji best¨atigt [TT91]. Die Bedeutung von Partikel-Partikel-Kollisionen kann, wie bereits in Kapitel 2 ausgef¨ uhrt, durch den Vergleich der Partikelrelaxationszeit τp mit der Kollisionszeit τC abgesch¨atzt werden. In einer dichten Gas-Feststoff-Str¨omung beispielsweise ist die Kollisionszeit deutlich kleiner als die Partikelrelaxationszeit. Hier wird die Partikelbewegung in entscheidendem Maße durch die interpartikul¨aren Kollisionen bestimmt. Die Modelle zur Simulation von Partikel-Partikel-Kollisionen unterscheiden sich in den Verfahren zur Detektion einer Kollision sowie in der Auswahl eines Kollisionspartners. Allen Modellen ist jedoch gemeinsam, dass die Ermittlung der Bewegungszust¨ande nach dem 30

3.2 Beschreibung der Feststoffphase Stoß aus der L¨osung der Impuls- und Drehimpuls-Erhaltungsgleichungen hervorgeht. Unter der Annahme, die Partikelphase bestehe aus starren Kugeln, werden dabei Inelastizit¨atsund Reibungseffekte u ¨ber den Restitutionskoeffizienten und den Gleitreibungsbeiwert in die Betrachtung miteinbezogen (Details siehe unter [Som96]). Zur Detektion einer Partikelkollision stehen sowohl deterministische als auch stochastische Modelle zur Verf¨ ugung. Bei einer deterministischen Kollisionsdetektion, wie es beispielsweise beim Modell von Tanaka & Tsuji der Fall ist, werden alle kombinatorisch denkbaren Partikelpaarungen in jedem Zeitschritt u uft [TT91]. Dies f¨ uhrt zu einem enormen Rechenaufwand, der ¨berpr¨ proportional zum Quadrat der Anzahl der Partikel ist. Daher ist diese Methode f¨ ur numerische Simulationen praxisrelevanter Zweiphasenstr¨omungen ungeeignet [CST98]. Um dennoch eine effiziente Berechnung von Gas-Feststoff-Str¨omungen unter Ber¨ ucksichtigung von Partikelkollisionen zu erm¨oglichen, wurden stochastische Kollisionsmodelle entwickelt. Bei diesen wird die Wahrscheinlichkeit f¨ ur das Auftreten einer Kollision abgesch¨atzt und die Kollisionsentscheidung anhand von Zufallszahlen getroffen. Ein Beispiel eines sehr genauen und effizienten, stochastischen Kollisionsmodells stellt die modifizierte Nanbu-Methode dar [Nan80, IN87]. Bei diesem Modell wird der Kollisionspartner zuf¨allig aus den N∆Ω Partikeln ausgew¨ahlt, die sich zu einem bestimmten Zeitpunkt t in der gleichen Zelle ∆Ω befinden. Die Wahrscheinlichkeit Pij f¨ ur das Auftreten einer Kollision zwischen dem i-ten und dem j-ten Partikel ist durch Pij =

νij ∆t N∆Ω

(3.61)

gegeben, wobei νij die Kollisionsfrequenz zwischen dem i-ten und j-ten Partikel darstellt und sich wie folgt berechnen l¨asst: νij = np π d2p |vij | .

(3.62)

Dabei bezeichnet np die lokale Anzahldichte der realen Feststoffpartikel. Die Wahrscheinlichkeit Pij entspricht der Aufenthaltswahrscheinlichkeit des j-ten Partikels im Kollisionszylinder, der vom i-ten und j-ten Partikel gebildet wird. Ein Kollisionszylinder, wie er schematisch in Abbildung 3.2 dargestellt ist, schließt alle geometrischen Orte ein, an denen innerhalb des Zeitschrittes ∆t eine Kollision zwischen den Partikeln i und j stattfinden kann. Hierzu muss sich der Schwerpunkt des j-ten Partikels im Inneren des Kollisionszylinders befinden, dessen Querschnittsfl¨ache sich bei monodispersen Partikeln zu πdp 2 ergibt. Die L¨ange des Zylinders entspricht dem Produkt aus der Relativgeschwindigkeit zwischen beiden Partikeln und der Zeitspanne ∆t. Der Kollisionspartner j wird mit Hilfe einer gleichverteilten Zufallszahl ψ ∈ [0, 1] ausgew¨ahlt, j = ψN∆Ω + 1,

(3.63)

wobei ψN∆Ω der ganzzahlige Anteil des Produkts ψN∆Ω ist. Das i-te Partikel kollidiert dann mit dem j-ten Partikel, falls folgende Bedingung erf¨ ullt ist: ψ>

j − Pij . N∆Ω

(3.64) 31

3 Formulierung der Bilanzgleichungen Hierbei sollte N∆Ω Pij < 1 gelten. Der Vorteil der modifizierten Nanbu-Methode liegt darin, dass der Rechenzeitbedarf direkt proportional zur Anzahl der Partikel in der betrachteten Zelle ist (O(N∆Ω)). Weiterhin besitzt der im Zufallsprozeß ausgew¨ahlte Kollisionspartner die Eigenschaften eines der simulierten Partikel innerhalb der Zelle ∆Ω. Dadurch wird sichergestellt, dass die Verteilungsfunktion der Partikeleigenschaften exakt repr¨asentiert wird. Dies verhindert im Gegensatz zu Kollisionsverfahren, die einen Kollisionspartner mit mittleren Partikeleigenschaften f¨ ur alle simulierten Partikel ausw¨ahlen ([OP93]), einen systematischen Modellfehler, der unter bestimmten Umst¨anden gravierende Auswirkungen auf die Str¨omungseigenschaften haben kann [Kan03]. Aus diesen Gr¨ unden kommt im Rahmen dieser Arbeit die modifizierte Nanbu-Methode zum Einsatz. Nach der Auswahl eines Kollisionspartners z und einer positiven Kollisionsentscheidung (Bedingung (3.64) ist erf¨ ullt) wird der Stoßpunkt j vij zwischen beiden Partikel bestimmt. Hierzu wird y das globale Koordinatensystem in das Koordinatensystem des Kollisionszylinders gem¨aß Abi bildung 3.2 transformiert. Anhand zweier gleichverteilter Zufallszahlen Γ1,2 ∈ [0, 1] erfolgt dann x die Erzeugung des Stoßnormaleneinheitsvektors auf einer Halbkugeloberfl¨ache des i-ten Parti¨ kels. Uber eine Koordinatentransformation wird der Stoßnormaleneinheitsvektor zur¨ uck in das Abbildung 3.2: Kollisionszylinder. globale Koordinatensystem transformiert. Aus der Position des Stoßnormaleneinheitsvektors ergibt sich schließlich der gesuchte Stoßpunkt zwischen den beiden Kollisionspartnern [Bab89].

3.3

Randbedingungen

In den vorangehenden Abschnitten 3.1 und 3.2 wurden die Bewegungsgleichungen f¨ ur beide Phasen einer Gas-Feststoff-Str¨omung vorgestellt. Die Formulierung der Randbedingungen vervollst¨andigt dieses Gleichungssystem und erm¨oglicht die Berechnung von konkreten Str¨omungsproblemen. Im Folgenden sollen jeweils f¨ ur das Zwei-Kontinua- und das Einzelpartikel-Verfahren die Randbedingungen an den Grenzen des Berechnungsgebietes f¨ ur die Gas- und Feststoffphase beschrieben werden.

3.3.1

Randbedingungen f¨ ur das Zwei-Kontinua-Verfahren

Randbedingungen f¨ ur die Gasphase An Eintrittsr¨andern wird die so genannte Dirichlet-Randbedingung verwendet, d.h. Werte aller Str¨omungsgr¨oßen der Gasphase werden fest vorgegeben. Am Austrittsrand werden die Gradienten in Fl¨achennormalenrichtung f¨ ur alle Variablen gleich Null gesetzt. An Symmetrier¨andern kann kein gerichteter Fluss einer Gr¨oße stattfinden, woraus folgt, dass 32

3.3 Randbedingungen ∂ug /∂n = 0 ist. Die Vorgabe eines Gradienten einer Gr¨oße am Geometrierand, wie es an Austritts- und an den Symmetrier¨andern erfolgt ist, wird als Neumann-Randbedingung bezeichnet. An festen W¨anden gilt die Haftbedingung f¨ ur die Geschwindigkeit der Gasphase, ug = 0. Die Turbulenzgr¨oßen an einer Wand betragen: ε = k = 0. Da diese Bedingung wegen der Notwendigkeit der Aufl¨osung des Wandnahbereiches nicht praktikabel ist, wird Folgendes gefordert:     ∂k ∂ε = 0, = 0, (3.65) ∂n W ∂n W wobei n die Koordinatenrichtung normal zur Wand ist. Um die turbulenten Gr¨oßen im wandnahen Bereich zu simulieren, kommen spezielle Wandfunktionen zum Einsatz. Hierzu wird der wandnahe Bereich in zwei Schichten aufgeteilt. In unmittelbarer Wandn¨ahe ist die Str¨omung haupts¨achlich durch die viskosen Spannungen gepr¨agt, man spricht in diesem Bereich von einer viskosen Unterschicht. Durch ¨ einen Ubergangsbereich, wo sowohl viskose Spannungen als auch Reynolds-Spannungen gleichermaßen bedeutsam sind, gelangt man in den Bereich der turbulenten Grenzschicht. Hier wird der Impulstransport von den Reynolds-Spannungen dominiert. Um den G¨ ultigkeitsbereich der beiden Schichten anzugeben, wird ein dimensionsloser Wandabstand Y + eingef¨ uhrt: Y+ =

ρg uτ y . µg

(3.66)

(3.66) enth¨alt die Wandschubspannungsgeschwindigkeit uτ , die sich aus der Wandschubspannung ergibt:
τW uτ = . (3.67) ρg Der Bereich der viskosen Unterschicht liegt zwischen Y + = 0 und Y + = 5. Der Bereich der turbulenten Grenzschicht erstreckt sich u ¨ber 50 ≤ Y + ≤ 5000, wobei die obere Grenze des G¨ ultigkeitsbereiches von der Reynolds-Zahl der Kernstr¨omung abh¨angig ist. Der ¨ Ubergangsbereich kann somit durch 5 ≤ Y + ≤ 50 beschrieben werden. Aufgrund der sehr geringen Ausdehnung der viskosen Unterschicht wird diese im Allgemeinen nicht durch das Rechengitter aufgel¨ost. Der wandn¨achste Punkt des numerischen ¨ Gitters liegt somit im Ubergangsbereich der turbulenten Grenzschicht und die zugeh¨orige Wandfunktion f¨ ur die Berechnung der mittleren Gasgeschwindigkeit kann wie folgt angegeben werden: |ug | =



uτ log Y + + 5.2, κ

(3.68)

wobei κ die von K´arm´an-Konstante mit κ = 0.4 darstellt. (3.68) wird als logarithmisches Wandgesetz bezeichnet. Um nun die Turbulenzgr¨oßen am wandn¨achsten Punkt zu berechnen, wird angenommen, dass sich die Produktion und die Dissipation der Turbulenz im 33

3 Formulierung der Bilanzgleichungen Gleichgewicht befinden. Unter dieser Annahme ergibt sich f¨ ur die Dissipationsrate der folgende Ausdruck [LS74]: ε=

u3τ . κy

(3.69)

Die turbulente kinetische Energie k l¨asst sich dann aus der L¨osung von (3.13) bestimmen. Randbedingungen f¨ ur die Feststoffphase Die Formulierung der Randbedingungen f¨ ur die Feststoffphase an Eintritts- und Austrittsr¨andern des Berechnungsgebietes erfolgt analog zur Gasphase. Als schwierig erweist sich die Behandlung der Randbedingungen an festen, undurchl¨assigen W¨anden. Die Problematik soll anhand der Formulierung der Randbedingung f¨ ur die Geschwindigkeit der Feststoffphase kurz skizziert werden. Die auf die Wand treffenden Partikel werden von der Wand selbst und von anderen Partikeln in Wandn¨ahe beeinflusst. Drei verschiedene Szenarien sind bei der Vorgabe der Randbedingungen f¨ ur die Geschwindigkeit an der Wand m¨oglich: • die Haftbedingung [CP95, HS97], wenn das Wandrauhigkeitsmaß gr¨oßer als der Partikeldurchmesser ist, die Partikel eine sehr geringe Tr¨agheit besitzen und somit eine vernachl¨assigbare Relativgeschwindigkeit zur Gasphase haben. In diesem Fall gilt: vp = 0.

(3.70)

Fast immer ist der Partikeldurchmesser gr¨oßer als die L¨angenskala der Wandrauigkeit, so dass Partikel an der Wand abgleiten k¨onnen. Auch bei einem so genannten Haftstoß, bei dem das Partikel schon w¨ahrend der Kompressionsphase des Stoßes auf der Wand abzurollen beginnt, ist die Haftbedingung nicht erf¨ ullt, weil es aufgrund eines Abrollens zu einer Relativbewegung zwischen der Wandfl¨ache und dem Partikel kommt, • die kinetische Gleichgewichtsbedingung [CSL86], d.h. die Partikel k¨onnen auf der Wand frei abgleiten und es tritt keine Differenz zwischen der mittleren Geschwindigkeit unmittelbar an der Wand und im Mittelpunkt der wandn¨achsten Zelle. Es gilt:   ∂up vp = 0, = 0, (3.71) ∂n W wobei die Gr¨oßen vp und up die Normal- und Tangentialgeschwindigkeit der Partikel an der Wand darstellen. Diese Bedingung ist nur dann sinnvoll, wenn der Einfluss der Wandreibung vernachl¨assigt werden kann,

34

3.3 Randbedingungen • die Schlupfbedingung, d.h. die Partikel werden an der Wand teilweise abgebremst. Es wird angenommen, dass die Schlupfgeschwindigkeit proportional zum lokalen Geschwindigkeitsgradienten ist [Soo67]:   ∂up up = lp , (3.72) ∂n W wobei lp ein L¨angenmaß kennzeichnet, f¨ ur den der mittlere Partikelabstand eingesetzt wird [DLS93]. (3.72) ist aus der kinetischen Theorie verd¨ unnter Gase abgeleitet worden und kann nur unter bestimmten Einschr¨ankungen auf makroskopische Partikel angewendet werden [Ind01]. Die hier vorgestellten Randbedingungen f¨ ur die wandparallelen Geschwindigkeitskomponenten der dispersen Phase ber¨ ucksichtigen nicht die mechanischen Vorg¨ange bei PartikelWand-Kollisionen. Eine k¨ urzlich von Indenbirken vorgestellte Methode l¨ost diese Problematik, indem die Partikelphase in zwei Teilphasen aufgespalten wird, die unterschiedliche Bewegungsrichtungen (eine Teilphase bewegt sich auf die Wand zu, die andere von der Wand weg) besitzen k¨onnen [Ind01]. Die mechanischen Vorg¨ange der Kollision eines Einzelpartikels finden dabei ihre Ber¨ ucksichtigung. Dieses Modell liefert neue Randbedingungen, die einen wesentlich allgemein g¨ ultigen Charakter aufweisen, als es bei den vorhergehenden Modellen der Fall war. Die Aufspaltung und Zusammenf¨ uhrung der Teilphasen f¨ uhrt jedoch zu einem h¨oherem Speicher- und Rechenaufwand, insbesondere dann, wenn die Anzahl der Teilphasen groß wird. Die oben dargelegte, offensichtliche Problematik bei der Formulierung von Randbedingungen f¨ ur die Geschwindigkeit der dispersen Phase tritt ebenfalls bei der Vorgabe von Randbedingungen f¨ ur die granulare Temperatur Θ auf. Im Rahmen dieser Arbeit wird im Bereich von festen, undurchl¨assigen W¨anden das Einzelpartikel-Verfahren eingesetzt. Aufgrund seiner detaillierten Modellierung stellt die Formulierung der Randbedingungen f¨ ur die disperse Phase an festen Grenzen des Berechnungsgebietes kein Problem mehr dar. Im Folgenden soll der Vorteil des Einzelpartikel-Verfahrens diesbez¨ uglich kurz erl¨autert werden.

3.3.2

Randbedingungen f¨ ur das Einzelpartikel-Verfahren

Die Formulierung der Randbedingungen f¨ ur die Gasphase bleibt bei der Verwendung des Einzelpartikel-Verfahrens f¨ ur die Feststoffphase unver¨andert. Aufgrund der individuellen Betrachtung eines jeden einzelnen Partikels der Feststoffphase ist die Behandlung der Randbedingungen f¨ ur die disperse Phase g¨anzlich unterschiedlich. Am Str¨omungseintritt werden gem¨aß entsprechender Verteilungsfunktionen eine Str¨omungsgeschwindigkeit sowie eine Partikelverteilung angenommen. Am Str¨omungsaustritt verlassen die Partikel ungehindert das Str¨omungsgebiet und werden in der Simulation nicht weiter ber¨ ucksichtigt. An festen, undurchl¨assigen W¨anden kollidieren die Partikel mit der Berandung. Im Gegensatz zu der kollektivistischen Betrachtungsweise im Zwei-Kontinua-Verfahren erm¨oglicht die Einzelpartikelbetrachtung eine detaillierte Beschreibung von Partikel-Wand-Wechselwirkungen. 35

3 Formulierung der Bilanzgleichungen Hierf¨ ur werden Informationen u ¨ber die mechanischen Eigenschaften der Kombination Partikel/Wand sowie u ¨ber den Bewegungszustand der Partikel vor der Kollision verwendet (s. Abschnitt 3.2.2). Durch die Anwendung der mechanischen Stoßgesetze zur Berechnung einer Wandkollision entf¨allt damit die Notwendigkeit der Formulierung von Randbedingungen f¨ ur die Feststoffphase an festen R¨andern. Aus diesem Grund wird im Rahmen dieser Arbeit in Bereichen fester Str¨omungsbegrenzungen ausschließlich das EinzelpartikelVerfahren verwendet.

3.4

Numerisches Berechnungsverfahren

Die in den Abschnitten 3.1 und 3.2.1 vorgestellten Gleichungen zur Berechnung der Gasphase und der dispersen Phase im Rahmen des Zwei-Kontinua-Verfahrens stellen ein System partieller, nichtlinearer Differentialgleichungen mit einer intensiven Kopplung der Variablen dar. Da eine analytische L¨osung des Gleichungssystems nicht m¨oglich ist, m¨ ussen effiziente, numerische Verfahren eingesetzt werden. Diese sind in der Literatur ausf¨ uhrlich dokumentiert [Fle88, FP96, Wen96] und sollen daher in dieser Arbeit nur kurz skizziert werden. Auf die L¨osung der Bewegungsgleichungen f¨ ur die disperse Phase im Rahmen des Einzelpartikel-Verfahrens (s. Abschnitt 3.2.2) soll am Ende dieses Abschnittes eingegangen werden. Zur numerischen L¨osung des Differentialgleichungssystems wird dieses in einem Diskretisierungsschritt durch ein algebraisches Gleichungssystem approximiert. Die Diskretisierung der Differentialgleichungen erfolgt nach der Finite-Volumen-Methode. Hierzu wird das Str¨omungsfeld in Kontrollvolumina unterteilt, die als Bilanzr¨aume zur L¨osung der integralen Transportgleichungen dienen. Alle Interpolationen und Gradientenapproximationen sind von zweiter Ordnung. Das verwendete Rechengitter ist blockstrukturiert, die Zellen stellen ausschließlich Hexaederelemente dar. Die Druck-Geschwindigkeit-Kopplung wird u ¨ber den SIMPLE-Consistent Algorithmus [VR84], ein Druckkorrektur-Verfahren sichergestellt. Das resultierende nichtlineare und gekoppelte Gleichungssystem wird iterativ gel¨ost. Die zur Zeit effizientesten Methoden zur iterativen L¨osung solcher Gleichungssysteme stellen die Mehrgitter-Verfahren dar. Dabei werden neben dem eigentlichen Rechengitter (Feingitter), auf dem die L¨osung des Problems erfolgen soll, zus¨atzlich eine Sequenz von Grobgittern mit einer geringeren Anzahl von Knotenpunkten erzeugt. Das Problem wird nun schrittweise auf die gr¨oberen Gitter u ¨bertragen, um dann auf dem gr¨obsten Gitter gel¨ost zu werden. Anschließend wird aus der Grobgitter-L¨osung und der Feingitter-Approximation eine Korrektur ermittelt, mit der die L¨osung auf dem urspr¨ unglichen Rechengitter verbessert wird. Eine detaillierte Beschreibung des eingesetzten Mehrgitter-Verfahrens erfolgt in der Arbeit von St¨ uben und Trottenberg [ST84]. Die Berechnung einer Gas-Feststoff-Str¨omung mit dem Einzelpartikel-Verfahren f¨ ur die disperse Phase erfolgt iterativ in zwei Teilschritten. Im ersten Schritt werden die Gleichungen zur Berechnung der Gasphase unter Anwendung des oben erw¨ahnten MehrgitterVerfahrens f¨ ur einen Zeitschritt ∆t gel¨ost. In einem zweiten Schritt erfolgt die Simulation 36

3.4 Numerisches Berechnungsverfahren der Feststoffphase. Um den Rechen- und Speicheraufwand zu begrenzen, werden im Allgemeinen nicht alle Partikel durch das Str¨omungsfeld verfolgt, sondern nur eine repr¨asentative Anzahl. Ein Simulationspartikel entspricht dann einer bestimmten Anzahl realer Partikel mit denselben physikalischen Eigenschaften. Bei der Simulation der Feststoffphase werden Wechselwirkungen sowohl mit der Gasphase als auch mit anderen Partikeln ber¨ ucksichtigt. Dabei m¨ ussen Zeitschrittbeschr¨ankungen eingehalten werden, die eine Aufteilung des Zeitschrittes ∆t in Subzeitschritte ∆τ erforderlich machen [Som96]. Bei der Ausf¨ uhrung eines Subzeitschrittes kommt die, im Abschnitt 3.2.2 bereits erw¨ahnte, Splitting-Technik zum Einsatz [Bird76]. Es handelt sich hier um eine Entkoppelung der Partikelbewegung von den interpartikul¨aren Kollisionen. Alle simulierten Partikel werden zun¨achst innerhalb von ∆τ kollisionsfrei bewegt. Treten w¨ahrend dieses Subzeitschrittes Wechselwirkungen mit der Berandung auf, so werden diese unmittelbar berechnet. Erst am Ende von ∆τ erfolgt die stochastische Ausf¨ uhrung der Partikel-Partikel-Kollisionen. Die Berechnung der Partikeltrajektorien wird solange fortgef¨ uhrt, bis die Summe aller Subzeitschritte den Zeitschritt ∆t der Gasphase erreicht. Nun k¨onnen die mittleren Eigenschaften der Partikelphase in jeder Zelle des Berechnungsgebietes ermittelt werden. Aus diesen Informationen werden schließlich die Quellterme in den Gleichungen der Gasphase, nach dem in Abschnitt 3.2.2 dargestellten PSIC-Verfahren, berechnet. Mit der erneuten Berechnung der Gasphase ist der erste Iterationszyklus abgeschlossen. Das iterative Verfahren wird solange fortgef¨ uhrt, bis ein Abbruchkriterium erf¨ ullt ist. Damit sind die numerischen Berechnungsmethoden f¨ ur die Berechnung einer Gas-Feststoff-Str¨omung jeweils unter Anwendung des Zwei-Kontinua- und des Einzelpartikel-Verfahrens f¨ ur die disperse Phase erl¨autert. Im Rahmen dieser Arbeit werden beide Verfahren zur Beschreibung der Feststoffphase miteinander kombiniert. Der Ablauf des gekoppelten Simulationsverfahrens wird in Kapitel 6 ausf¨ uhrlich beschrieben.

37

3 Formulierung der Bilanzgleichungen

38

4

Analyse des Verhaltens der Feststoffphase

Wie im vorigen Kapitel erl¨autert, l¨asst sich das Bewegungsverhalten eines Molek¨ ulkollektivs in einem idealen Gas durch die Angabe einer Wahrscheinlichkeitsdichteverteilung f (t, x, v) beschreiben. Wirken auf das System keinerlei Ausgleichsprozesse infolge treibender Kr¨afte, so befindet sich das Molek¨ ulkollektiv im thermodynamischen Gleichgewicht. In diesem Zustand entspricht die Wahrscheinlichkeitsdichteverteilung der bekannten MaxwellVerteilung. Das Bewegungsverhalten eines Partikelkollektivs in einer turbulenten Gasstr¨omung ist trotz seiner Analogie zur Bewegung eines Molek¨ ulkollektivs erheblich komplexer. Der Gleichgewichtszustand entspricht nicht dem eines idealen Gases. Die Abweichung vom thermodynamischen Gleichgewicht resultiert aus der Existenz zus¨atzlicher Einflußfaktoren, die in der Natur einer Feststoffstr¨omung liegen (s. Abschnitt 3.2.1). In diesem Kapitel soll ausgehend von der Definition des Gleichgewichtszustandes f¨ ur ein Partikelkollektiv, ein Gleichgewichtskriterium vorgestellt werden, welches das Ausmaß der Abweichung einer Verteilungsfunktion von der Gleichgewichtsverteilung bewertet. Auf der Grundlage dieses Kriteriums sollen dann Richtlinien f¨ ur die Anwendbarkeit des Zwei-Kontinua- bzw. Einzelpartikel-Verfahrens zur Berechnung der Feststoffphase formuliert werden.

4.1

Das granulardynamische Gleichgewicht

Betrachtet man ein Partikelkollektiv in einem unendlich ausgedehnten Raum, so stellt sich nach unendlich langer Zeit f¨ ur die Wahrscheinlichkeitsdichteverteilung der Partikelgeschwindigkeiten f (v) ein Gleichgewichtszustand ein. Dieser ist dadurch charakterisiert, dass sich die Verteilungsfunktion nicht mehr ¨andert und somit keine Ausgleichsprozesse mehr stattfinden. F¨ ur diesen Zustand hat Kanther in seiner Arbeit den Begriff des granulardynamischen Gleichgewichts eingef¨ uhrt [Kan03]. Befindet sich ein Partikelkollektiv im granulardynamischen Gleichgewicht, so unterscheidet sich die zugeh¨orige Verteilungsfunktion der Partikelgeschwindigkeiten f (v) von der Maxwell-Verteilung durch zwei charakteristische Deformationen. Zum einen durch die Anisotropie, die aus dem Einfluss der Gravitationskraft resultiert, zum anderen durch eine strukturelle Abweichung in der Asymptotik infolge inelastischer Partikel-Partikel-Kollisionen. Bei der Anisotropie handelt es sich um eine achsensymmetrische Deformation der Maxwell-Verteilung, d.h. die Schwankungsgeschwindigkeiten der Partikel besitzen eine unterschiedliche Auspr¨agung in den Koordinatenrichtungen. Die strukturelle Abweichung in der Asymptotik einer Verteilungsfunktion spielt eine wichtige Rolle nur im Bereich kleiner Stokes-Zahlen und f¨ ur sehr geringe Restitutionskoeffizienten. Betrachtet man dreidimensionale, berandete Gas-Feststoff-Str¨omungen, so ist der Einfluss inelastischer Kollisionen auf die Geschwindigkeitsverteilung der Partikel vernachl¨assig39

4 Analyse des Verhaltens der Feststoffphase bar klein. Die Gr¨ unde hierf¨ ur liegen darin, dass die Kollisionen nahezu elastisch erfolgen und bei h¨oheren Partikel-Reynoldszahlen (Rep > 1) von der Anisotropie der Str¨omung dominiert werden. Da in technisch relevanten Gas-Feststoff-Str¨omungen deutlich h¨ohere Partikel-Reynoldszahlen zu erwarten sind, wird die Geschwindigkeitsverteilung maßgeblich von der Anisotropie beeinflusst. Damit l¨asst sich das granulardynamische Gleichgewicht durch eine angen¨aherte Gleichgewichtsverteilung fEq,p charakterisieren, die anisotrope Geschwindigkeitsfluktuationen bereits ber¨ ucksichtigt [Kan03]:   3  1 (vi − ui )2 √ fEq,p (v) = np exp − . (4.1) 2Θi 2πΘi i=1 Der Unterschied zur Maxwell-Verteilung (3.18) liegt in der Verwendung richtungsabh¨angiger granularer Temperaturen Θi . Bisher wurde Θ unter der Annahme isotroper Geschwindigkeitsfluktuationen als eine richtungsunabh¨angige Gr¨oße verwendet. Diese Annahme ist beim Vorliegen anisotroper Fluktuationen nicht mehr zul¨assig. Zur Charakterisierung der Abweichungen von (4.1) l¨asst sich eine Bewertungsfunktion formulieren, die in folgendem Abschnitt vorgestellt wird.

4.2

Erkennung granulardynamischer Gleichgewichtszust¨ ande

Um das Ausmaß der Abweichung einer beliebigen Verteilungsfunktion fp der Partikelgeschwindigkeiten von der Gleichgewichtsverteilung (4.1) bewerten zu k¨onnen, kann eine Bewertungsfunktion DEq,p formuliert werden, f¨ ur deren Herleitung die Eigenschaften beider Verteilungsfunktionen fp und fEq,p verglichen werden. Als Eigenschaften eignen sich statistische Momente, die sich f¨ ur beide Funktionen explizit bestimmen lassen. Aus den Differenzen dieser Gr¨oßen und unter Verwendung eines geeigneten Normierungsmaßstabes kann dann eine skalare Gr¨oße gefunden werden, welche eine Beurteilung der Verteilungsfunktion fp der Partikelgeschwindigkeiten erm¨oglicht. Bevor ein entsprechendes Kriterium zur Beurteilung von Verteilungsfunktionen gegen¨ uber dem granulardynamischen Gleichgewicht vorgestellt wird, sollen zun¨achst zwei Methoden skizziert werden, mit denen eine Bewertungsfunktion angegeben werden kann, welche die Abweichungen einer beliebigen Verteilungsfunktion fp der Partikelgeschwindigkeiten von der Maxwell-Verteilung fM erfasst. Hierzu wird jeweils ein r¨aumlich homogenes System von N Feststoffpartikeln betrachtet. Die erste Methode stellt eine exakte M¨oglichkeit zur Berechnung des Abstandes DEq einer Verteilungsfunktion fp zur lokalen MaxwellVerteilung fM dar. Hierbei wird DEq unter Verwendung der so genannten Sobolev-Norm bestimmt [TR97]:  DEq := fM − fp  = 40

R3

|fˆM (ω) − fˆp (ω)|2 dω (1 + |ω|2 )

 12 ,

(4.2)

4.2 Erkennung granulardynamischer Gleichgewichtszust¨ande wobei die Gr¨oßen fˆM (ω) bzw. fˆp (ω) Fourier-Transformierte im Frequenzraum ω ∈ R3 darstellen. (4.2) erfasst alle Deformationen einer Verteilungsfunktion gegen¨ uber einer MaxwellVerteilung. Im Folgenden soll die Gesamtheit der auf eine Verteilungsfunktion einwirkenden Deformationen als Typ-I-Deformationen bezeichnet werden. Ergibt die Ermittlung von (4.2) einen Wert gleich Null, so liegt mit Sicherheit eine Maxwell-Verteilung vor, anderenfalls existieren strukturelle Abweichungen von dieser. Der Nachteil dieser Methode ist, dass zur Auswertung von (4.2) ein Aufwand der Ordnung O(N 2) erforderlich ist, weshalb diese Methode f¨ ur die Praxis ungeeignet ist. In einer weiteren von Tiwari vorgestellten Methode wird f¨ ur die zu untersuchende Verteilungsfunktion fp eine Abweichung von der Maxwell-Verteilung durch den folgenden Zusammenhang ber¨ ucksichtigt [Tiw98a, Tiw98b]: fp (v) = fM (v) [1 + χ(v)] .

(4.3)

F¨ ur die Funktion χ(v) wird ein polynomialer Ansatz nach der Grad´schen 13 Momente Methode [Grad49] gemacht, die Kanther f¨ ur die Beschreibung von Partikelkollektiven um ein Moment vierter Ordnung erweitert hat [Kan03]: χ(v) = a + b · (v − u) + (v − u)t C (v − u) + d · (v − u) |v − u|2 + e · |v − u|4 . (4.4) Gleichung (4.4) enth¨alt 14 Unbekannte, von denen a und e skalare Gr¨oßen, b und d Vektoren und C einen symmetrischen Tensor charakterisieren. Zur Ermittlung der unbekannten Koeffizienten werden statistische Momente nullter bis vierter Ordnung der diskreten Verteilungsfunktion bestimmt: m(0) =

N 

gj ,

m(1) =

j=1

m(3) =

N  j=1

N 

gj vj ,

m(2) =

N 

j=1

gj vj |vj |2 ,

m(4) =

gj vj vjT ,

(4.5)

j=1 N 

gj |vj |4 .

j=1

Hierbei bezeichnet gj das Partikelgewicht. Mit Hilfe dieses Gewichtungsfaktors muss die Anzahl an simulierten Partikeln im Str¨omungsgebiet nicht mit der Anzahl real existierender Partikel u ¨bereinstimmen. Aus den statistischen Momenten lassen sich nun die ben¨otigten makroskopischen Gr¨oßen bestimmen. Die Anzahldichte np entspricht dabei dem Quotienten aus dem Moment nullter Ordnung und dem Volumen einer Zelle: np =

m(0) . ∆Ω

(4.6)

Die drei Komponenten der mittleren Partikelgeschwindigkeit k¨onnen unter Einbeziehung des Momentes erster Ordnung gebildet werden: u=

m(1) . m(0)

(4.7) 41

4 Analyse des Verhaltens der Feststoffphase Der symmetrische Spannungstensor T, die granulare Temperatur Θ und der W¨armeflußvektor q ergeben sich aus: T = m(2) − m(0) uuT ,

Θ=

3 1  Tii , 3m(0) i=1

(4.8)

  3 1 1  (2) q = m(3) − m(2) u + u m(0) |u|2 − mii . 2 2 i=1 F¨ ur die Bestimmung der Bewertungsfunktion DEq ist die Angabe zweier weiterer Gr¨oßen erforderlich, des spurfreien Spannungstensors τ und des zentralen Momentes vierter Ordnung γ: τ = T − m(0) ΘI,

(4.9)

γ = m(4) − 8 u · m(3) + 4 uT m(2) u + 2 m(2) I |u|2 − 3 |u|4 − 15m(0) Θ2 . Unter Verwendung von (4.8) und (4.9) sowie einer Normierungsfunktion f¨ ur den polynomialen Ansatz χ(v)  DEq ≈ χ :=

R3

fM (v) 2 χ (v) dv ρ

 12 (4.10)

ur die Distanz zwischen der erh¨alt man schließlich die Bewertungsfunktion DEq , als ein Maß f¨ Maxwell-Verteilung fM und einer beliebigen Verteilung fp der Partikelgeschwindigkeiten [Kan03]: DEq ≈ χ =



1 1 γ2 2 |q|2 τ 2E + + (0) m Θ 2 5 Θ 120 Θ2 1

 12 .

(4.11)

Hierbei bezeichnet ||τ ||E die Euklidische Norm des spurfreien Spannungstensors. Die in der eckigen Klammer von (4.11) enthaltenen Summanden beschreiben drei unterschiedliche Deformationsarten. Anisotrope Abweichungen von der Maxwell-Verteilung werden durch den Spannungstensor ||τ ||E erfasst, asymmetrische Deformationen in der Verteilung der kinetischen Energie ber¨ ucksichtigt der W¨armeflußvektor q. Den strukturellen Abweichungen in der Asymptotik einer Verteilungsfunktion fp wird schließlich durch das zentrale Moment vierter Ordnung γ Rechnung getragen. Diese drei Deformationsarten werden im Folgenden als Typ-II-Deformationen bezeichnet. Gleichung (4.11) stellt eine Bewertungsfunktion zur Beurteilung einer beliebigen Verteilungsfunktion fp der Partikelgeschwindigkeiten gegen¨ uber der Maxwell-Verteilung fM dar. Der Aufwand f¨ ur die Auswertung von (4.11) ist von der Ordnung O(N ). Um Abweichungen einer Verteilungsfunktion fp von der granulardynamischen Gleichgewichtsverteilung fEq,p detektieren zu k¨onnen, muss die Bewertungsfunktion DEq durch die Verwendung richtungsabh¨angiger granularer Temperaturen Θi modifiziert werden. Die neue Bewertungsfunktion 42

4.3 Granulare Transportkoeffizienten DEq,p nimmt dann die folgende Form an: DEq,p =

  12 1 2 2 2 1 1 Dτ + Dq + Dγ2 , (0) m 2 5 120

(4.12)

wobei die einzelnen Summanden Dτ , Dq und Dγ wiederum die Deformationsarten Anisotropie, Asymmetrie und Asymptotik bezeichnen. Eine ausf¨ uhrliche Darstellung dieser Terme findet sich in der Arbeit von Kanther [Kan03]. Mit (4.12) k¨onnen Abweichungen vom granulardynamischen Gleichgewicht zuverl¨assig erfasst werden, solange diese zum Typ-II-Deformationen geh¨oren. Da (4.12) in jeder Zelle des Berechnungsgebietes ermittelt wird, k¨onnen innerhalb der Rechengeometrie Partikelkollektive, die sich in einem granulardynamischen Nicht-Gleichgewichts-Zustand befinden, lokalisiert werden. Diese Bereiche werden als Nicht-Gleichgewichts-Zonen bezeichnet.

4.3

Granulare Transportkoeffizienten

Eine detaillierte Beschreibung der Verteilungsfunktion der Partikelgeschwindigkeiten ist von großer Bedeutung f¨ ur die korrekte Berechnung des Bewegungsverhaltens eines Partikelkollektivs. In diesem Zusammenhang sind drei Gr¨oßen zu nennen, die bei der Str¨omungsanalyse eine zentrale Position einnehmen. Es handelt sich dabei um die mittlere Kollisionsfrequenz ν, die mittlere Fluktuationsgeschwindigkeit c und die mittlere freie Wegl¨ange
, f¨ ur die folgender Zusammenhang gilt:
=

c . ν

(4.13)

W¨ahrend die mittlere Kollisionsfrequenz ν die mittlere Anzahl der St¨oße eines Partikels in einer Zeiteinheit angibt und somit maßgeblich f¨ ur die pr¨azise Beschreibung des Kollisionsverhaltens der Partikel ist, sind die Gr¨oßen c und
wichtig f¨ ur die Herleitung der Transportkoeffizienten im Rahmen des Zwei-Kontinua-Verfahrens. Um die Abh¨angigkeit der drei Gr¨oßen von der Gestalt der Geschwindigkeitsverteilung zu untersuchen, hat Kanther entsprechende Parameterstudien durchgef¨ uhrt [Kan03]. Diese haben gezeigt, dass bereits geringe Abweichungen der Verteilungsfunktion fp von der Gleichgewichtsverteilung fEq,p zu erheblichen Abweichungen der Kollisionsfrequenz ν gegen¨ uber der Kollisionsfrequenz ν M einer energiegleichen Maxwell-Verteilung f¨ uhren k¨onnen. In diesem Zusammenhang ist auch die Wahl eines geeigneten Kollisionsmodells diskutiert worden. Es konnte aufgezeigt werden, dass die zur Detektion der Partikelkollisionen verwendete Nanbu-Methode die mittlere Kollisionsfrequenz ν eines Partikelkollektivs richtig wiedergibt. F¨ ur die Untersuchung der Abh¨angigkeit der granularen Transportkoeffizienten von der Form der Verteilungsfunktion fp wurde die Viskosit¨at µp eines Partikelkollektivs ausgew¨ahlt. Die Untersuchungen haben gezeigt, dass die Ber¨ ucksichtigung der Anisotropie im granulardynamischen Gleichgewicht zu einer anisotropen Viskosit¨at µp f¨ uhrt. Die Unterschiede in den gerichteten Komponenten der Viskosit¨at nehmen stark zu, 43

4 Analyse des Verhaltens der Feststoffphase wenn sich das Partikelkollektiv vom granulardynamischen Gleichgewicht entfernt. Dabei k¨onnen die Unterschiede zum Teil Gr¨oßenordnungen betragen. Dar¨ uber hinaus weichen die absoluten Gr¨oßen der Komponenten der Viskosit¨at betr¨achtlich von der Viskosit¨at µM eines Partikelkollektivs mit einer energiegleichen Maxwell-Verteilung ab. Diese Ergebnisse sind auch auf andere granulare Transportkoeffizienten u ¨bertragbar.

4.4

Schlussfolgerung

In diesem Kapitel wurde eine Bewertungsfunktion DEq,p vorgestellt, welche alle Abweichungen einer Verteilungsfunktion der Partikelgeschwindigkeiten fp von der granulardynamischen Gleichgewichtsverteilung fEq,p zuverl¨assig erfasst, solange diese eine anisotrope, asymmetrische oder asymptotische Deformation der Verteilungsfunktion zur Folge haben. Damit l¨asst sich eine beliebige Verteilungsfunktion eines Partikelkollektivs hinsichtlich ihres Abstandes vom Gleichgewichtszustand bewerten. Befindet sich ein Partikelkollektiv nicht in einem granulardynamischen Gleichgewicht, so hat dies entscheidende Konsequenzen f¨ ur die Ermittlung der Transportkoeffizienten, die im Rahmen des Zwei-Kontinua-Verfahrens ben¨otigt werden. Bereits bei der Annahme einer granulardynamischen Gleichgewichtsverteilung weisen die Transportkoeffizienten einen anisotropen Charakter auf. Dieser wird verst¨arkt, sobald sich das Partikelkollektiv vom Gleichgewicht entfernt. In diesem Fall werden nicht nur die Unterschiede in den gerichteten Transportkoeffizienten wesentlich gr¨oßer, sondern sie weichen auch betr¨achtlich von den Transportkoeffizienten ab, die eine energiegleiche Maxwell-Verteilung aufweist. Das Zwei-Kontinua-Verfahren, welches sich aus der deterministischen L¨osung der Boltzmann-Gleichung ergibt (vgl. Abschnitt 3.2.1), setzt an jedem Ort des Berechnungsgebietes und zu jeder Zeit die Maxwell-Verteilung sowie geringe Abweichungen von dieser voraus. Wegen der Annahme einer isotropen Geschwindigkeitsverteilung, besteht bei der Herleitung der Transportkoeffizienten ein prinzipieller Modellfehler. Dieser wird jedoch vernachl¨assigbar klein, wenn Partikelkollektive mit h¨oheren Volumenanteilen untersucht werden. Hier sind die Abweichungen aufgrund der Anisotropie meistens nur geringf¨ ugig. Bei gravierenden Deformationen der Geschwindigkeitsverteilung jedoch, kann das Bewegungsverhalten eines Partikelkollektivs im Rahmen des Zwei-Kontinua-Verfahrens nicht richtig wiedergegeben werden. Solche Deformationen der Verteilungsfunktion treten insbesondere in komplexen Geometrien auf, beispielsweise vor Str¨omungshindernissen, wo die Annahme einer isotropen Geschwindigkeitsverteilung zur Berechnung der Transportkoeffizienten auf keinen Fall gerechtfertigt ist. In diesen Nicht-Gleichgewichts-Zonen f¨ uhrt die Berechnung einer Partikelstr¨omung mit dem Zwei-Kontinua-Verfahren zu fehlerhaften Ergebnissen. ¨ Die oben formulierten Uberlegungen legen den Schluss nahe, dass die Verwendung des Zwei-Kontinua-Verfahrens auf Bereiche innerhalb der Str¨omungsgeometrie begrenzt werden sollte, in denen sich das Partikelkollektiv in einem Zustand nahe dem granulardynamischen Gleichgewicht befindet. Außerhalb dieser Bereiche, in denen die Verteilungsfunktion der Partikelgeschwindigkeiten eine starke Deformation erf¨ahrt, ist die Verwendung des stochastischen Einzelpartikel-Verfahrens notwendig. Hier werden keine Einschr¨ankungen 44

4.4 Schlussfolgerung bez¨ uglich der Verteilungsfunktion gemacht, d.h. alle Abweichungen vom granulardynamischen Gleichgewicht fließen in die Auswertung der mittleren Partikeleigenschaften ein und werden insbesondere durch die Nanbu-Methode ber¨ ucksichtigt. Damit stellt sich aber auch die Frage nach welchem Kriterium eine Aufteilung der Str¨omungsgeometrie erfolgen soll, in der eine unterschiedliche numerische Behandlung der Berechnung der Feststoffphase zum Einsatz kommt. Dieser Frage widmet sich das folgende Kapitel, in dem die Herleitung eines Gebietszerlegungskriteriums vorgestellt wird.

45

4 Analyse des Verhaltens der Feststoffphase

46

5

Gebietszerlegungskriterium

Bei der Berechnung des Bewegungsverhaltens eines Partikelkollektivs mit dem Zwei-Kontinua-Verfahren werden f¨ ur die Verteilungsfunktion der Partikelgeschwindigkeiten im granulardynamischen Gleichgewicht eine Maxwell-Verteilung vorausgesetzt sowie bestimmte Abweichungen von dieser zugelassen. Weitere Abweichungen von diesem Zustandsraum f¨ uhren zu einer fehlerhaften Berechnung der Transportkoeffizienten und damit zu einer unrealistischen Vorhersage des Partikelverhaltens. Wie im vorigen Kapitel erl¨autert, muss in Bereichen, in denen die Deformationen der Geschwindigkeitsverteilung stark anisotrope Transportkoeffizienten zur Folge haben, das stochastiHindernis sche Einzelpartikel-Verfahren eingesetzt werden. Diese Bereiche werden im Folgenden als stochastische Zonen bezeichnet. Im restlichen Teil der Geometrie kann das Zwei-Kontinua-Verfahren verstochastische Zone wendet werden. Diese Bereiche werden als deterministische Zonen bezeichnet. Die Entscheidung f¨ ur die Verwendung eines der beiden numedeterministische Zone rischen Berechnungsverfahren kann mit Hilfe eines Gebietszerlegungskriteriums erfolgen. Dieses Abbildung 5.1: Gebietszerlegung. sollte die zeitliche Evolution der Eigenschaften des Partikelkollektivs ber¨ ucksichtigen und damit eine adaptive Berechnungsweise erm¨oglichen, ohne den G¨ ultigkeitsbereich des Zwei-Kontinua-Verfahrens zu verletzen. In Abbildung 5.1 ist am Beispiel der Umstr¨omung eines Hindernisses schematisch die Aufteilung der Str¨omungsgeometrie in zwei Bereiche dargestellt, in denen die Berechnung des Partikelkollektivs auf numerisch unterschiedlichem Wege erfolgt. In diesem Kapitel werden zwei Bedingungen vorgestellt, die eine Aufteilung der Geometrie in Gebiete erm¨oglichen, in denen das jeweils geeignete Berechnungsverfahren zum Einsatz kommt. Diese m¨ ussen dann durch logische Verkn¨ upfungen zu einem Gebietszerlegungskriterium zusammengefasst werden.

5.1

Gu ¨ ltigkeitsbereich der Kontinuumsannahme

Der Berechnung einer Gas-Feststoff-Str¨omung mit dem Zwei-Kontinua-Verfahren liegt die Annahme zugrunde, dass sich sowohl die Gas- als auch die Feststoffphase als zwei kontinuierliche Fluide verhalten, sich gegenseitig durchdringen und miteinander in Wechselwirkung stehen. Die Grundvoraussetzung f¨ ur eine Kontinuumsbeschreibung beider Phasen ist, dass f¨ ur einen Mittelwert eines beliebigen Kontrollvolumens des durchstr¨omten Gebietes die mittleren, lokalen Zustandsgr¨oßen jeweils in geeigneter und statistisch sinnvoller Weise definiert werden k¨onnen. Daraus resultiert die Forderung, dass die Anzahl von Einzelelemen47

5 Gebietszerlegungskriterium ten beider Phasen f¨ ur jedes betrachtete Kontrollvolumen im Str¨omungsgebiet hinreichend groß sein muss. F¨ ur die Feststoffphase bedeutet dies insbesondere, dass der mittlere Partikelabstand deutlich geringer sein sollte als die Abmessungen des gew¨ahlten numerischen ¨ Gitters. Zur Uberpr¨ ufung der Kontinuumsannahme ist es u ¨blich, statt des mittleren Partikelabstandes die mittlere freie Wegl¨ange
zu verwenden, f¨ ur die folgende Absch¨atzung verwendet werden kann [EPA99]:
=

dp . 6 φp

(5.1)

Es handelt sich hier um die mittlere Wegstrecke, die zwischen zwei aufeinander folgenden Kollisionen eines Partikels mit einem anderen Partikel oder einer Berandung zur¨ uckgelegt wird. Mit der mittleren freien Wegl¨ange
und einer charakteristischen Abmessung lchar eines Kontrollvolumens l¨asst sich nun eine dimensionslose Kennzahl, die Knudsen-Zahl Kn definieren, welche als Bedingung f¨ ur die G¨ ultigkeit der Kontinuumshypothese verwendet werden kann: Kn =


≤ Kncrit , lchar

(5.2)

√ wobei lchar z.B. aus dem Volumen einer Gitterzelle berechnet werden kann: lchar = 3 V∆Ω . Schaaf&Chambr` e teilen die Knudsen-Zahl in drei Bereiche ein (s. Abbildung 5.2), in denen eine unterschiedliche mathematische Beschreibung eines Gases erfolgt [SC61]. So

10−8

10−4

10−6

Kontinuums-Bereich

10−2

¨ Ubergangs-Bereich

1

102

104

Kn

freie Molek¨ ulstr¨omung

Abbildung 5.2: Bereiche der Knudsen-Zahl. ist die G¨ ultigkeit der Kontinuumsbetrachtung bis zu Kncrit = 10−3 gew¨ahrleistet. F¨ ur −1 Kn ≥ 10 ist f¨ ur eine statistische Mittelung die Anzahl der Molek¨ ule in einem Kontrollvolumen nicht ausreichend groß, weshalb hier eine Einzelbetrachtung der Molek¨ ule erfol¨ gen sollte (freie Molek¨ ulstr¨omung). Zwischen den beiden Bereichen, im Ubergangs-Bereich (10−3 ≤ Kn ≤ 10−1 ), muss im Einzelfall entschieden werden, welche Berechnungsmethode zum Einsatz kommt. W¨ahrend die Gasphase die Bedingung Kn ≤ 10−3 in jedem Fall erf¨ ullt, ist sie f¨ ur die in dieser Arbeit betrachteten Partikelkollektive nicht realisierbar. Dies gilt insbesondere dann, wenn die charakteristische L¨ange lchar sehr klein wird. Es hat sich jedoch gezeigt, dass viele praxisrelevante Mehrphasenstr¨omungen in sehr guter N¨aherung mit dem Zwei-Kontinua-Verfahren beschrieben werden k¨onnen. Der Grund daf¨ ur liegt darin, dass 48

5.2 KTGF-Deformationsanteil DKT GF zumeist station¨are oder zumindest quasistation¨are Str¨omungen betrachtet werden. In solchen F¨allen k¨onnen nach Goedicke die Ensemble-Mittelwerte durch Langzeitmittelwerte ersetzt werden [Goe92]. Folglich l¨asst sich die Kontinuumsannahme durch die Verwendung hinreichend langer Mittelungszeitr¨aume rechtfertigen. Um dennoch bei der Berechnung der Bewegung eines Partikelkollektivs mit dem Zwei-Kontinua-Verfahren eine Mindestanzahl an Partikeln in einer Zelle des numerischen Gitters zu gew¨ahrleisten, wird f¨ ur die Simulation der Feststoffphase die Knudsen-Zahl wie folgt festgelegt: Kncrit = 10−1 . Damit kann als Bedingung f¨ ur die Anwendbarkeit des Zwei-Kontinua-Verfahrens die folgende Grenze f¨ ur den Feststoffvolumenanteil angegeben werden: φp >

dp . 0.6 lchar

(5.3)

Sinkt der Feststoffvolumenanteil in einem Kontrollvolumen unter den angegebenen Wert ab, so kann eine repr¨asentative Anzahl an Partikeln f¨ ur eine instation¨are Kontinuumsbetrachtung nicht gew¨ahrleistet werden. Somit muss die Beschreibung des Bewegungsverhaltens des Partikelkollektivs in dem betrachteten Kontrollvolumen mit dem EinzelpartikelVerfahren approximiert werden.

5.2

KTGF-Deformationsanteil DKT GF

Die in dem vorigen Kapitel vorgestellte Gleichgewichtsverteilung fEq,p der Partikelgeschwindigkeiten (4.1) beschreibt den Bewegungszustand eines Partikelkollektivs, das sich in einem granulardynamischen Gleichgewicht befindet. Abweichungen von diesem Zustand werden mit Hilfe der Bewertungsfunktion DEq,p (4.12) detektiert. Dabei werden alle St¨orungen der Verteilungsfunktion fp zuverl¨assig erfasst, die als eine anisotrope, asymmetrische oder asymptotische Deformation der Geschwindigkeitsverteilung aufgefasst werden k¨onnen (Typ-II-Deformationen). Damit kann die Bewertungsfunktion DEq,p als ein Unterscheidungskriterium zwischen einer Gleichgewichts- und einer Nicht-Gleichgewichts-Verteilung der Partikelgeschwindigkeiten verwendet werden. Beim Zwei-Kontinua-Verfahren werden im granulardynamischen Gleichgewicht bestimmte Abweichungen einer Verteilungsfunktion von der Maxwell-Verteilung ber¨ ucksichtigt. Um den Abweichungsgrad genauer zu charakterisieren, wird ein, urspr¨ unglich von Chapman & Cowling vorgestellter [CC70] und von Gidaspow auf Partikelkollektive u ¨bertragener Ansatz zur Beschreibung dichter Gase verwendet [Gid94], der eine Approximation der Verteilungsfunktion der Partikelgeschwindigkeiten darstellt: f = f (0) + f (1) + f (2) + .... .

(5.4)

f (0) bezeichnet die erste Approximation zur Boltzmann-Gleichung und stellt die bekannte Maxwell-Verteilung dar. f = f (0) ist somit eine erste Approximation der Geschwindigkeitsverteilung und gilt zun¨achst f¨ ur r¨aumlich homogene Partikelkollektive mit vollkommen elastischen Partikel-Partikel-Kollisionen. Abweichungen von diesem Zustand werden 49

5 Gebietszerlegungskriterium durch weitere Glieder in (5.4) ber¨ ucksichtigt. So bezeichnen f (0) + f (1) die zweite und (0) (1) (2) f + f + f die dritte Approximation zur Boltzmann-Gleichung. Zur Ermittlung von f (1) in (5.4) wird folgender Ansatz verwendet: f (1) = f (0) Φ(1) ,

(5.5)

wobei Φ(1) ein skalares Abweichungspolynom darstellt, welches bestimmte Deformationen der Maxwell-Verteilung erfasst. Unter Verwendung der Boltzmann-Gleichung (3.16) kann f¨ ur Φ(1) der folgende Zusammenhang gefunden werden [Gid94]: Φ(1) = −

1 2 A · ∇lnΘ − B : ∇uTp , np np

(5.6)

wobei A eine vektorielle und B eine tensorielle Gr¨oße darstellen. Beide sind Funktionen der Partikelgeschwindigkeit v. F¨ ur A wird folgender Ansatz verwendet: A = A(ξ) ξ,

(5.7)

wobei ξ eine dimensionslose Variable mit  ξ=

1 2Θ

 12 v

(5.8)

darstellt und die Funktion A(ξ) eine konvergente Reihe der Form A(ξ) =

∞  p=0

(p)

ap S 3 (ξ 2 )

(5.9)

2

ist. Im Rahmen dieser Arbeit wird die unendliche Reihe in (5.9) durch die zweite Approximation von A(ξ) ersetzt: A(ξ) ≈ A(ξ)(2) =

2  p=0

(p)

2 a(2) p S 3 (ξ ).

(5.10)

2

(2)

Dabei erfolgt die Ermittlung der Koeffizienten ap aus der Auswertung folgender Gleichung: 2 

a(2) p apq = αq , (q = 1, 2).

(5.11)

p=0

Hierbei sind α1 = −15/4 und α2 = 0. Die fehlenden Koeffizienten apq k¨onnen [CC70] (p) entnommen werden. In (5.9) bezeichnet S 3 (ξ 2 ) das so genannte Sonin-Polynom, das sich 2 aus einer Reihendarstellung ergibt:   p  2 k 3 1 (p)

+p . (5.12) −ξ S 3 ξ2 = 2 2 k!(p − k)! p−k k=0 50

5.2 KTGF-Deformationsanteil DKT GF Die Auswertung von (5.11) sowie (5.12) f¨ uhrt schließlich zu folgender Beziehung f¨ ur A(ξ): A(ξ) ≈ −

15 130 − 59ξ 2 + 2ξ 4 . 1 1408 d2p (π Θ) 2

(5.13) ◦

Zur Ermittlung der tensoriellen Gr¨oße B wird das Produkt aus einem spurfreien Tensor ξξ und einer konvergenten Reihe B(ξ) gebildet: ◦

B =ξξ B(ξ).

(5.14)

◦

Der Tensor ξξ ist definiert als: ◦

ξξ= ξξ −

1 (ξx ξx + ξy ξy + ξz ξz ) I, 3

(5.15)

wobei ξξ das dyadische Produkt der beiden ξ-Vektoren darstellt. Die Berechnung von B(ξ) erfolgt in analoger Weise zu A(ξ) aus der Auswertung von B(ξ)(2) =

2 

(p−1)

(ξ 2 )

(5.16)

b(2) p bpq = βq , (q = 1, 2),

(5.17)

p=0

b(2) p S5 2

sowie 2  p=0

ur B(ξ) ergibt sich danach der folgende Zusamwobei β1 = 5/2 und β2 = 0 betragen. F¨ menhang : B(ξ) ≈ −

5 −247 + 12ξ 2 . 3232 d2 (π Θ) 12 p

(5.18)

Mit (5.6) ist es nun m¨oglich zwei unterschiedliche Ursachen f¨ ur die Deformationen einer Verteilungsfunktion anzugeben. Der erste Summand in (5.6) ber¨ ucksichtigt Abweichungen von der Maxwell-Verteilung aufgrund eines Gradienten der granularen Temperatur. Der zweite Summand erm¨oglicht die Erfassung aller Abweichungen, die aus dem Gradienten der Partikelgeschwindigkeit resultieren. Damit ist die zweite Approximation der Verteilungsfunktion der Partikelgeschwindigkeiten in der Lage Verteilungsfunktionen von Partikelkollektiven richtig zu beschreiben, solange deren Form durch die beiden genannten Gradienten deformiert worden ist. Um weitere Abweichungen von der Gleichgewichtsverteilung zu ber¨ ucksichtigen, welche beispielsweise aufgrund h¨oherer Ableitungen entstehen, m¨ ussen h¨ohere Approximationen der Verteilungsfunktion in die Berechnung einbezogen werden. Da f¨ ur die Herleitung der Transportgleichungen des Zwei-Kontinua-Verfahrens die 51

5 Gebietszerlegungskriterium

zweite Approximation der Verteilungsfunktion f = f (0) 1 + Φ(1) verwendet wird (vgl. Abschnitt 3.2.1), kann im Weiteren auf eine Auswertung entsprechender Abweichungspolynome h¨oherer Ordnung verzichtet werden. Mit der Angabe des Abweichungspolynoms (5.6) lassen sich nun folgende Aussagen formulieren. Kommt es zu einer Deformation der Verteilungsfunktion aufgrund von Gradienten der granularen Temperatur oder der Partikelgeschwindigkeit, so werden diese im Rahmen des Zwei-Kontinua-Verfahrens erfasst. Treten jedoch zus¨atzliche Abweichungen der Verteilungsfunktion von der Gleichgewichtsverteilung auf, die ihren Ursprung nicht in den genannten St¨orungen haben, so bleiben diese innerhalb des genannten Verfahrens unber¨ ucksichtigt. In diesem Fall w¨ urde die Berechnung des Partikelverhaltens zu systematisch falschen Ergebnissen f¨ uhren. Das Abweichungspolynom Φ(1) und die im letzten Kapitel vorgestellte Bewertungsfunktion DEq,p sollen nun dazu dienen eine Bedingung zu formulieren, die eine Aufteilung der Str¨omungsgeometrie in stochastische und deterministische Zonen erm¨oglicht, in denen jeweils das Einzelpartikel- und das Zwei-Kontinua-Verfahren zum Einsatz kommen. Wie bereits oben ausgef¨ uhrt, ber¨ ucksichtigt das Abweichungspolynom Φ(1) Deformationen aus Gradienten der Geschwindigkeit und der granularer Temperatur. Die Bewertungsfunktion DEq,p dagegen ist in der Lage alle Deformationsarten aufgrund von Anisotropie, Asymmetrie und Asymptotik zu erfassen. Diese schließen ebenfalls die vom Abweichungspolynom Φ(1) beschriebenen Deformationen ein. Somit stellt die Bewertungsfunktion DEq,p ein zuverl¨assiges Maß dar, um den aus den erw¨ahnten Deformationen resultierenden Abweichungsgrad einer Verteilungsfunktion vom granulardynamischen Gleichgewicht zu bewerten. Die Idee einer Gebietszerlegung besteht nun darin, das Abweichungspolynom Φ(1) mit dem Maßstab DEq,p zu bewerten. Um die beiden Gr¨oßen Φ(1) und DEq,p miteinander vergleichen zu k¨onnen, muss (5.6) mit einer geeigneten Normierungsfunktion verrechnet werden. Die hier verwendete Norm entspricht der f¨ ur die Herleitung der Bewertungsfunktion beschriebenen Darstellung [Tiw98a, Tiw98b]:  χ :=

R3

fM (v) 2 χ (v) dv ρ

 12 ,

(5.19)

in welcher f¨ ur die Funktion χ(v) zuvor ein polynomialer Ansatz nach der Grad´schen 13 Momenten Methode gemacht wurde. Wird χ(v) durch das Abweichungspolynom Φ(1) ersetzt, so ergibt sich der folgende Ausdruck:  DKT GF =

R3

f (0) (v) (1) 2 Φ dv np

 12 .

(5.20)

Setzt man nun (5.6) in (5.20) ein, so l¨asst sich f¨ ur DKT GF nach einigen Umformungen die nachstehende Formulierung finden: DKT GF = 52

0.763 1

Θ np dp 2 π 2

|∇Θ| +

0.184 2

np dp (Θ π)

1

1 2

(2A1 + A2 + 3A3 ) 2 ,

(5.21)

5.2 KTGF-Deformationsanteil DKT GF wobei folgende Abk¨ urzungen gelten:  2  2  2 ∂up ∂vp ∂wp A1 = + + , ∂x ∂y ∂z  A2 =  A3 =

∂up ∂vp − ∂x ∂y ∂up ∂vp + ∂y ∂x



2 +



2 +

∂up ∂wp − ∂x ∂z ∂up ∂wp + ∂z ∂x

(5.22) 

2 +



2 +

∂vp ∂wp − ∂y ∂z ∂vp ∂wp + ∂z ∂y

2 ,

(5.23)

.

(5.24)

2

Gleichung (5.21) stellt den KTGF-Deformationsanteil dar, der alle Abweichungen von der Maxwell-Verteilung erfasst, solange diese ihren Ursprung in den Gradienten der granularen Temperatur und der Partikelgeschwindigkeit haben. Dieser Typ von Deformationen wird im Folgenden als Typ-III-Deformationen bezeichnet. Zus¨atzliche Deformationen der Verteilungsfunktion werden von DKT GF nicht ber¨ ucksichtigt. Es ist nun m¨oglich den KTGF-Deformationsanteil DKT GF (5.21) direkt mit der Bewertungsfunktion DEq,p (4.12) zu vergleichen. In der Abbildung 5.3 sind beide Gr¨oßen in einem Balkendiagramm dargestellt. Ausgehend vom granulardynamischen Gleichgewicht DEq,p , DKT GF Typ-I-Deformationen DEq,p (Typ-II-Deformat.)

Asymptotik Abweichungen von der KTGF Asymmetrie

DKT GF (Typ-III-Deformat.)

∇Θ ∇Θ

Anisotropie

∇vT EinzelpartikelVerfahren

∇vT Zwei-KontinuaVerfahren

granulardynamisches Gleichgewicht

Abbildung 5.3: Vergleich der Bewertungsfunktion DEq,p mit dem KTGF-Deformationsanteil DKT GF . zeigt der linke Balken die Gesamtheit aller Deformationen, denen eine Verteilungsfunktion der Partikelgeschwindigkeiten ausgesetzt werden kann. In der Abbildung 5.3 wird dieser 53

5 Gebietszerlegungskriterium Zustand durch die gestrichelte, horizontale Linie mit dem Index: Typ-I-Deformationen gekennzeichnet. Abweichungen vom granulardynamischen Gleichgewicht, welche von der Bewertungsfunktion DEq,p erfasst werden, sind durch die horizontale Linie mit dem Index: DEq,p (Typ-II-Deformationen) gekennzeichnet. Hierin enthalten sind die Anisotropie, die Asymmetrie und die Asymptotik. Die anisotropen St¨orungen der Verteilungsfunktion k¨onnen zu einem Teil aus den Gradienten der Partikelgeschwindigkeit (∇vT ) resultieren. Der restliche Teil, hier durch die untere, schraffierte Fl¨ache dargestellt, enth¨alt weitere St¨orungen, wie sie beispielsweise vor Str¨omungshindernissen entstehen. Sie werucksichtigt. Anders stellt sich die Situation f¨ ur den in der Bewertungsfunktion DEq,p ber¨ den KTGF-Deformationsanteil DKT GF dar, der durch den rechten Balken charakterisiert wird. Ausschließlich St¨orungen, die aus dem Gradienten der Partikelgeschwindigkeit re¨ sultieren, tragen zu einer anisotropen Deformation der Verteilungsfunktion bei. Ahnlich verh¨alt es sich beim Vergleich der asymmetrischen Deformationen einer Verteilungsfunktion. So ber¨ ucksichtigt der KTGF-Deformationsanteil DKT GF nur Abweichungen aufgrund des Gradienten der granularen Temperatur ∇Θ. Weitere St¨orungen der Verteilungsfunktion, die zur einer asymmetrischen Deformation f¨ uhren, werden im Gegensatz zur DEq,p , hier durch die schraffierte, mittlere Fl¨ache dargestellt, nicht ber¨ ucksichtigt. Eine asymptotische Deformation einer Verteilungsfunktion wird schließlich nur von der Bewertungsfunktion DEq,p detektiert, im KTGF-Deformationsanteil findet sie keine Ber¨ ucksichtigung. Die Tatsache, dass der KTGF-Deformationsanteil ausschließlich Typ-III-Deformationen erfasst, f¨ uhrt dazu, dass das Zwei-Kontinua-Verfahren eine gute Approximation des Bewegungsverhaltens eines Partikelkollektivs liefert, solange die vorliegenden Deformationen zu diesem Deformationstyp geh¨oren. In der Abbildung (5.3) ist dieser Zustand durch die gestrichelte, horizontale Linie mit dem Index: DKT GF (Typ-III-Deformationen) gekennzeichnet. Tragen weitere Abweichungen zur Deformation der Verteilungsfunktion bei, so werden diese vom DKT GF nicht detektiert und f¨ uhren somit bei einer Berechnung der Feststoffphase mit dem Zwei-Kontinua-Verfahren zur fehlerhaften Ergebnissen. In diesen F¨allen ist der Einsatz des Einzelpartikel-Verfahrens zwingend erforderlich. Zus¨atzliche Abweichungen, die nicht zum Typ-III-Deformationen geh¨oren, sind in der Abbildung 5.3 als Abweichungen von der KTGF gekennzeichnet und werden zu einem großen Teil von der Bewertungsfunktion DEq,p detektiert. Die verbleibenden Deformationen, in der Abbildung (5.3) durch die obere, schraffierte Fl¨ache dargestellt, geh¨oren zum Typ-I-Deformationen. Deren Erfassung kann mit Hilfe der in Kapitel 4 beschriebenen Sobolev-Norm erfolgen. Damit kann eine weitere Bedingung formuliert werden, deren G¨ ultigkeit den zwingenden Einsatz des Einzelpartikel-Verfahrens fordert: DEq,p > DKT GF .

(5.25)

Im Falle, dass DEq,p = DKT GF kann das Zwei-Kontinua-Verfahren angewendet werden.

54

5.3 Signifikanz-Niveau-Beschr¨ankung

5.3

Signifikanz-Niveau-Beschr¨ ankung

Bei der Simulation des Bewegungsverhaltens eines Partikelkollektivs werden die beiden Gr¨oßen DEq,p und DKT GF nach jedem Zeitschritt ∆t und in jeder Zelle der stochastischen Zone ermittelt. Folglich unterliegt die Berechnung beider Gr¨oßen stochastischen Fluktuationen. W¨ahrend diese f¨ ur die Ermittlung des KTGF-Deformationsanteils DKT GF eine untergeordnete Rolle spielen, hier wird ausschließlich auf die in den Zellenmittelpunkten vorliegenden Gradienten zur¨ uckgegriffen, k¨onnen sich bei der Berechnung der Bewertungsfunktion DEq,p sehr hohe Unsicherheiten ergeben. Diese sind umso gr¨oßer, je kleiner die Anzahl der zur Auswertung ben¨otigten Partikel ist. Folglich kann die Ermittlung der im vorigen Abschnitt vorgestellten Bedingung (5.25) zu einer falschen Entscheidung hinsichtlich der Wahl des numerischen Berechnungsverfahrens f¨ uhren. Um dies zu vermeiden, wird die Bewertungsfunktion DEq,p mit einer Funktion verglichen, die sich aus einem numerischen Experiment ergibt. Hierzu werden r¨aumlich homogene Partikelkollektive im granulardyna0,9 0.9

DEq,p,hom

0.6 0,6

Mittelwert D Eq,p,hom Standardabweichung σ(DEq,p,hom )

0,3

0.3

0

0

200 200

400 400

np

600 600

800 800

1000 1000

Abbildung 5.4: Signifikanz-Niveau der Bewertungsfunktion DEq,p . mischen Gleichgewicht erzeugt, die jeweils die gleiche Anzahl an Partikeln besitzen. Deren Simulation ergibt im station¨aren Endzustand eine Bewertungsfunktion DEq,p,hom . Aus der Vielzahl dieser Simulationen kann schließlich f¨ ur eine definierte Anzahl von Partikeln ein Mittelwert D Eq,p,hom und eine Standardabweichung σ(DEq,p,hom ) ermittelt werden. Dieses Experiment wird f¨ ur verschiedene Partikelzahlen fortgef¨ uhrt. In der Abbildung 5.4 sind die beiden Parameter in Abh¨angigkeit der Partikelanzahl np dargestellt. Es sind bereits die an die Simulationsdaten angepassten Ausgleichsgeraden eingezeichnet. Es ergibt sich der folgende Zusammenhang zwischen der mittleren Bewertungsfunktion DEq,p,hom bzw. deren Standardabweichung σ(DEq,p,hom ) und der Partikelanzahl np : DEq,p,hom = 2.795 n−0.493 , p σ(DEq,p,hom ) =

0.835 n−0.525 . p

(5.26) (5.27) 55

5 Gebietszerlegungskriterium Bei sehr großen Partikelzahlen streben beide Gr¨oßen asymptotisch einem sehr kleinen Wert entgegen. Damit l¨asst sich nun zusammen mit (5.25) die folgende Bedingung f¨ ur die Anwendbarkeit eines der beiden numerischen Berechnungsverfahren in einer beliebigen Zelle ∆Ω der Rechengeometrie formulieren: DEq,p > DKT GF + 1.96 σ(DEq,p,hom ),

(5.28)

wobei hier ein 95%-Konfidenzintervall angenommen wird. Ist (5.28) erf¨ ullt, so liegt eine stochastische Zone vor und das Einzelpartikel-Verfahren kommt in jedem Fall zum Einsatz.

5.4

Gebietszerlegungskriterium

Mit Hilfe der in den Abschnitten 5.1 und 5.3 vorgestellten Bedingungen (5.3) und (5.25) sowie unter Ber¨ ucksichtigung von (5.28) kann das folgende Gebietszerlegungskriterium formuliert werden. F¨ ur den Einsatz des Einzelpartikel-Verfahrens gilt: ∗ DEq,p > DKT GF

oder

Kn > Kncrit ,

(5.29)

∗ wobei DKT ur den GF = DKT GF + 1.96 σ(DEq,p,hom ). Es liegt eine stochastische Zone vor. F¨ Einsatz des Zwei-Kontinua-Verfahrens gilt: ∗ DEq,p ≤ DKT GF

und

Kn ≤ Kncrit .

(5.30)

In diesem Fall liegt eine deterministische Zone vor.

5.5

Schlussfolgerung

In diesem Kapitel wurde ein Kriterium (5.29) und (5.30) vorgestellt, mit dessen Hilfe eine Entscheidung bez¨ uglich der Anwendung eines der beiden numerischen Berechnungsverfahren getroffen wird. Das Kriterium basiert auf zwei Bedingungen. W¨ahrend die erste (5.3) die Kontinuumsannahme und damit die G¨ ultigkeit des Zwei-Kontinua-Verfahrens u uft, ¨berpr¨ nutzt die zweite Bedingung (5.25) die entscheidenden Annahmen aus, die f¨ ur die Herleitung des Zwei-Kontinua-Verfahrens aus der Boltzmann-Gleichung zugrunde gelegt worden sind. In diesem Zusammenhang ist der so genannte KTGF-Deformationsanteil entwickelt worden, mit dessen Hilfe die Erfassung bestimmter Deformationen und ein direkter Vergleich mit der Bewertungsfunktion DEq,p m¨oglich sind. Die Ber¨ ucksichtigung von statistischen Schwankungen, denen die Auswertung der Bedingung (5.25) ausgesetzt ist, erfolgt durch die Einbeziehung einer Signifikanzfunktion. Die Auswertung des Kriteriums in jeder Zelle des Berechnungsgebietes erm¨oglicht in jedem Zeitschritt der Simulation eine Aufteilung der Str¨omungsgeometrie in stochastische und deterministische Zonen. In den deterministischen Zonen ist eine aggregierte Beschreibung des Einzelpartikelverhaltens eines Partikelkollektivs mit dem Zwei-Kontinua-Verfahren zul¨assig. In den stochastischen Zonen dagegen sind detaillierte Informationen u ¨ber die Verteilung der Partikelgeschwindigkeiten erforderlich, 56

5.5 Schlussfolgerung um das Bewegungsverhalten der Feststoffphase zuverl¨assig zu beschreiben. Hier ist der Einsatz des Einzelpartikel-Verfahrens zwingend erforderlich. Durch die Zuordnung eines numerischen Berechnungsverfahrens zu bestimmten Bereichen der Str¨omungsgeometrie ist es nun m¨oglich komplexe Str¨omungsprobleme unabh¨angig von dem Verhalten der sich darin befindlichen Partikel zu untersuchen. Als Voraussetzung daf¨ ur m¨ ussen allerdings Techniken entwickelt werden, die eine adaptive Berechnung der gesamten Str¨omungsgeometrie mit einer Kombination beider Einzelverfahren erm¨oglichen. Im Mittelpunkt des n¨achsten Kapitels steht daher die Konstruktion eines neuen, gekoppelten Verfahrens, welches die Berechnungsstrategien der beiden Einzelverfahren integriert und damit eine effiziente Simulation eines erweiterten Anwendungsbereichs erm¨oglicht.

57

5 Gebietszerlegungskriterium

58

6

Konstruktion eines gekoppelten, adaptiven Simulationsverfahrens

Im vorherigen Kapitel wurde ein Gebietszerlegungskriterium vorgestellt, das eine Aufteilung der Berechnungsgeometrie in zwei Bereiche erm¨oglicht, in denen eine unterschiedliche numerische Behandlung der Feststoffphase stattfinden soll. Die Ermittlung des Bewegungsverhaltens eines Partikelkollektivs erfolgt in beiden Geometriebereichen unter Verwendung der in Kapitel 3 vorgestellten Modellans¨atze f¨ ur das Einzelpartikel- bzw. das Zwei-Kontinua-Verfahren. Da es sich bei einer turbulenten Gas-Feststoff-Str¨omung prinzipiell um einen instation¨aren Vorgang handelt, muss nach jedem Iterationsschritt im Verlaufe der Simulation eine erneute Auswertung des in dem letzten Kapitel gefundenen Gebietszerlegungskriteriums (5.29) und (5.30) erfolgen. Dies impliziert, dass die r¨aumliche Ausdehnung der Bereiche in denen ein bestimmtes numerisches Verfahren verwendet wird, ¨ variieren kann und somit eine kontinuierliche Uberwachung der Gebietsaufteilung nicht nur numerisch sinnvoll, sondern auch aus Gr¨ unden der Modellkonsistenz erforderlich ist. In diesem Kapitel soll zun¨achst das Konzept der adaptiven Zerlegung des Berechnungsgebietes vorgestellt werden. Im Weiteren werden Techniken behandelt, die einen Austausch der Informationen zwischen den numerischen Verfahren erm¨oglichen. Anschließend wird die L¨osungstrategie des gekoppelten Simulationsverfahrens erl¨autert.

6.1

Konzept der adaptiven Verfahrenskopplung

Die Anwendung des im vorigen Kapitel vorgestellten Gebietszerlegungskriteriums erm¨oglicht eine eindeutige Zuordnung einzelner Geometriebereiche zu einem der beiden numerischen Berechnungsverfahren. Eine Aufteilung der Str¨omungsgeometrie in zwei Zonen, d.h. in eine stochastische Zone (S) in der das Einzelpartikel-Verfahren und eine zweite, deterministische Zone (D) in der das Zwei-Kontinua-Verfahren verwendet wird (s. Abbildung 6.1), w¨ urde je¨ doch eine automatische Uberwachung der Gebietsaufteilung verhindern. Der Grund daf¨ ur stochastische Zone (Einzelpartikel-Verfahren) deterministische Zone (Zwei-Kontinua-Verfahren)

S

∗ DEq,p > DKT GF ∨ Kn > Kncrit

D

∗ DEq,p ≤ DKT GF ∧ Kn ≤ Kncrit

Abbildung 6.1: Gebietszerlegung in zwei Zonen. liegt darin, dass die Auswertung des Gebietszerlegungskriteriums ausschließlich in Berei59

6 Konstruktion eines adaptiven Simulationsverfahrens chen m¨oglich ist, in denen das Einzelpartikel-Verfahren eingesetzt wird. Eine r¨aumliche Ausdehnung der stochastischen Zone (S) w¨are damit nicht feststellbar. Eine L¨osung dieses Problems erm¨oglicht die Definition einer Randzone (R). Diese besitzt die Breite von mindestens einer Zelle und ist dadurch charakterisiert, dass hier das Einzelpartikel-Verfahren zum Einsatz kommt, obwohl die Auswertung des Gebietszerlegungskriteriums eine deterministische Zone anzeigt (s. Abbildung 6.2). Der Vorteil dieser Methode ist, dass die Randzone S

∗ DEq,p > DKT GF ∨ Kn > Kncrit

Einzelpartikel-Verfahren R ∗ DEq,p ≤ DKT GF ∧ Kn ≤ Kncrit

Zwei-Kontinua-Verfahren

D

Abbildung 6.2: Gebietszerlegung in drei Zonen. nach jedem Zeitschritt r¨aumlich neu angepasst werden kann. Damit ist sichergestellt, dass die Aufteilung der Str¨omungsgeometrie in stochastische (S) und deterministische Zonen (D) nur von der Auswertung des Gebietszerlegungskriteriums abh¨angig ist und w¨ahrend der Simulation vollst¨andig adaptiv erfolgt. Ergibt die Auswertung des Gebietszerlegungs∗ kriteriums in einer Randzone: DEq,p > DKT GF oder Kn > Kncrit , so wird aus der Randzone (R) eine stochastische Zone (S). Umgekehrt kann aus der Randzone eine deterministische Zone (D) werden, wenn in allen benachbarten stochastischen Zonen (S) die Auswertung ∗ des Gebietszerlegungskriteriums DEq,p ≤ DKT GF und Kn ≤ Kncrit ergibt. Hierbei ist zu beachten, dass eine stochastische Zone (S) niemals an eine deterministische Zone (D) grenzen kann, d.h. sie muss stets von Randzonen (R) umgeben sein. In der Abbildung 6.3 ist die r¨aumliche Evolution der Gebietsaufteilung w¨ahrend einer Simulation f¨ ur die stochastische Zone (S) schematisch dargestellt. Die Auswertung des Gebietszerlegungskriteriums erfolgt nach jedem Zeitschritt und in jeder Zelle, die zur Rand- bzw. zur stochastischen Zone geh¨ort. Muss eine R-Zelle, in der Abbildung 6.3 links oben mit  R gekennzeichnet, zu einer S-Zelle werden, so expandiert die stochastische Zone und alle benachbarten D-Zellen, ¨ mit  D gekennzeichnet, der neuen S-Zelle werden zu R-Zellen (Expansion). Ahnliches gilt f¨ ur die Kontraktion der stochastischen Zone; soll aus einer S-Zelle, in der Abbildung 6.3 rechts unten mit  S gekennzeichnet, eine D-Zelle werden, so wechselt der Status dieser Zelle aufgrund der Nachbarschaft einer S-Zelle zun¨achst zu einer R-Zelle. Alle benachbarten R-Zellen, die ihrerseits keine S-Zelle als Nachbar haben, werden zu D-Zellen (Kontraktion). Nach jeder Ver¨anderung der Zellenzugeh¨origkeit muss sichergestellt sein, dass die stochastische Zone vollst¨andig von Randzellen umgeben ist. Auf diese Weise k¨onnen even¨ tuelle Anderungen des Bewegungsverhaltens eines Partikelkollektivs erkannt und in der oben beschriebenen Prozedur der Gebietszerlegung ber¨ ucksichtigt werden. ¨ Den bisherigen Uberlegungen ist gemeinsam, dass sie f¨ ur eine Erweiterung bzw. eine Verringerung der stochastischen Zone stets eine r¨aumliche Verschiebung der Randzone 60

6.1 Konzept der adaptiven Verfahrenskopplung Expansion D R

D D

D R

D

D

D R

R

D

R

S

R

D

R

S

S

D

R

S

S

R

D R

D

D

D R

D

D

D R

R

D

Kontraktion D D R

D

D R

R

D

D R

R

D

D R

S

R

R

R

S

R

R

S

S

R

D D R

D

D R

R

D

D R

R

D

Abbildung 6.3: Expansion und Kontraktion der stochastischen Zone (S). voraussetzen. Im Folgenden soll gekl¨art werden, ob eine stochastische Zone auch innerhalb eines Bereiches entstehen kann, der vollst¨andig zur deterministischen Zone geh¨ort. Die Ausbildung einer stochastischen Zone verlangt die Erf¨ ullung mindestens einer der bei∗ ultigkeit der ersten Bedinden Bedingungen: DEq,p > DKT GF oder Kn > Kncrit . Die G¨ ∗ gung (DEq,p > DKT orung des Bewegungsverhaltens des Partikelkollektivs GF ) setzt eine St¨ voraus, die vom KTGF-Deformationsanteil nicht erfasst wird. Eine solche St¨orung kann jedoch bei einer station¨aren Simulation nur u ¨ber Randbedingungen aufgepr¨agt werden. Diese k¨onnen im Verlaufe einer Simulation durch die Pr¨asenz von festen Berandungen entstehen. Da in dem hier verwendeten Algorithmus sichergestellt wird, dass s¨amtliche wandn¨achsten Zellen zur stochastischen Zone geh¨oren, k¨onnen etwaige St¨orungen durch die erste Bedingung erfasst werden. Somit kann innerhalb einer deterministischen Zone der ∗ Fall DEq,p > DKT GF gar nicht auftreten. Die Ausbildung einer stochastischen Zone aufgrund der G¨ ultigkeit der ersten Bedingung kann daher ausgeschlossen werden. In diesem Zusammenhang sei darauf hingewiesen, dass die Erfassung von stochastischen Inseln bereits zu Beginn einer Simulation erfolgen muss. Diese haben die Funktion einer Keimzelle und k¨onnen im Verlaufe der Simulation ihre endg¨ ultige Gestalt annehmen. Auf diese Weise k¨onnen Nicht-Gleichgewichts-Zonen, die beispielsweise eine Folge der sich gegenseitig durchdringenden Partikelkollektive sind, erfasst werden. Aus der Erf¨ ullung der zweiten Bedingung Kn > Kncrit folgt, dass der Volumenanteil der Feststoffphase in der betrachteten Zelle so gering ist, dass die Kontinuumsannahme f¨ ur die Anwendung des Zwei-Kontinua-Verfahrens nicht mehr erf¨ ullt ist. In diesem Fall muss die betreffende Zelle in eine stochastische Zelle umgewandelt werden. Bisher erfolgt die Auswertung des Gebietszerlegungskriteriums nur in den stochastischen- und den Randzo∗ ¨ nen. Der Grund daf¨ ur ist, dass die Uberpr¨ ufung der ersten Bedingung (DEq,p > DKT GF ) die Kenntnis u ur die ¨ber die Verteilungsfunktion der Partikelgeschwindigkeiten erfordert. F¨ Auswertung der zweiten Bedingung Kn > Kncrit wird lediglich der Volumenanteil der Fest61

6 Konstruktion eines adaptiven Simulationsverfahrens D

D

D

D

D

D

D

R

D

D

D

D

D

S

D

R

S

R

D

D

D

D

D

D

D

R

D

Abbildung 6.4: Ausbildung einer S-Zelle innerhalb einer deterministischen Zone (D) aufgrund der G¨ultigkeit der Bedingung Kn > Kncrit . stoffphase ben¨otigt, der im gesamten Rechengebiet zur Verf¨ ugung steht. Somit kann die Auswertung der zweiten Bedingung auch in den zur deterministischen Zone geh¨orenden DZellen erfolgen. In der Abbildung 6.4 ist der Vorgang der Ausbildung einer stochastischen Zelle innerhalb einer deterministischen Zone schematisch dargestellt. Ein Ablaufdiagramm der gekoppelten Simulation einer Gas-Feststoff-Str¨omung wird im Abschnitt 6.3 beschrieben.

6.2

Informationsaustausch an den Zonengrenzen

Die Ermittlung der Str¨omungseigenschaften eines Feststoffkollektivs mit einem hybriden Simulationsverfahren setzt voraus, dass sowohl das Einzelpartikel-Verfahren als auch das Zwei-Kontinua-Verfahren in einem Bereich, in dem dessen Anwendung zul¨assig ist, vergleichbare Ergebnisse liefern. Hierzu ist ein Vergleich beider Verfahren notwendig, um eventuelle Unterschiede in der Vorhersage des Partikelverhaltens zu lokalisieren. Ein weiteres wichtiges Element in der Entwicklung des neuen Verfahrens stellt der Informationsaustausch zwischen den beiden Einzelverfahren an der Zonengrenze dar. Es m¨ ussen ¨ entsprechende Rand- und Ubergangsbedingungen formuliert werden, die eine konservative Behandlung der Bewegungsgleichungen im gesamten Rechengebiet gew¨ahrleisten. Die oben genannten Aspekte sollen nun in den folgenden Abschnitten genauer behandelt werden.

6.2.1

Korrektur der Widerstandskraft

Die translatorische Bewegung eines Einzelpartikels wird in entscheidendem Maße von der Widerstandskraft FW gepr¨agt: FW = cd

ρg πdp 2 |vrel |vrel , 2 4

vrel = ug − vp .

(6.1)

Die Ermittlung des Widerstandsbeiwertes cd erfolgt dabei in Abh¨angigkeit der ReynoldsZahl der Partikel Rep . F¨ ur Rep < 0.25 spricht man vom Stokes-Bereich; auf die Partikel wirken vorwiegend Reibungskr¨afte. F¨ ur den Widerstandsbeiwert und die Widerstandskraft

62

6.2 Informationsaustausch an den Zonengrenzen gelten in diesem Bereich folgende Beziehungen: cd =

24 Rep

⇒

FW = 3π dp µ vrel .

(6.2)

F¨ ur Rep > 1000 sind es die Tr¨agheitskr¨afte, die einen entscheidenden Einfluss auf die Partikel aus¨ uben. In diesem so genannten Newton-Bereich liegt eine turbulente Umstr¨omung der Partikel vor und die Widerstandskraft ist proportional zum Quadrat der Relativgeschwindigkeit: cd = 0.44,

FW = 0.44

ρg π d2p vrel |vrel |. 8

(6.3)

¨ Im Ubergangsbereich 0.25 ≤ Rep ≤ 1000 sind sowohl Reibungs- als auch Tr¨agheitskr¨afte gleichermaßen relevant. In Kapitel 3 wurde zur Berechnung der Widerstandskraft eine empirische Korrelation angegeben (3.40), welche die in den drei Bereichen vorherrschenden Kraftverh¨altnisse ber¨ ucksichtigt. Wie (6.2) und (6.3) zeigen, ist die Widerstandskraft von der Relativgeschwindigkeit abh¨angig. Dabei ist vrel keine konstante Gr¨oße, sondern durch eine Verteilungsfunktion bestimmt. Schwankungen der Relativgeschwindigkeit werden somit auf die Widerstandskraft u ¨bertragen. Im Stokes-Bereich ergibt sich ein proportionaler Zusammenhang zwischen den Varianzen der Widerstandskraft σF2 W und der Relativgeschwindigkeit σv2 rel [Kan03]: σF2 W = (3π dp µ)2 σv2 rel .

(6.4)

Hierin ist σF2 W unabh¨angig von der mittleren Relativgeschwindigkeit vrel = ug − up . Im ¨ Ubergangssowie im Newton-Bereich ergibt sich dagegen ein nichtlinearer Zusammenhang zwischen σF2 W und σv2 rel . Zudem ist σF2 W zunehmend von der mittleren Relativgeschwindigkeit vrel abh¨angig. F¨ ur den Newton-Bereich ergibt sich zwischen σF2 W und σv2 rel bzw. vrel n¨aherungsweise folgende Abh¨angigkeit: 2 

ρg σF2 W ≈ 2 0.44 π d2p σv2 rel σv2 rel + 2 |vrel |2 . (6.5) 8 (6.5) bedeutet, dass aufgrund der zus¨atzlichen Abh¨angigkeit der Varianz der Widerstandskraft von vrel die Schwankungen der Relativgeschwindigkeit anisotrop auf die Schwankungen der Widerstandskraft u ¨bertragen werden. Daraus resultiert eine wichtige Konsequenz f¨ ur die Berechnung der Widerstandskraft. W¨ahrend das Zwei-Kontinua-Verfahren die Widerstandskraft aus der mittleren Relativgeschwindigkeit vrel ermittelt, FW (Rep ) = cd (Rep )

ρg πdp 2 |vrel |vrel , 2 4

Rep =

vrel ρg dp , µg

(6.6)

verwendet das Einzelpartikel-Verfahren die vollst¨andige Verteilung der Relativgeschwindigur ein Partikelkollektiv folgende mittlere Widerstandskraft ergibt: keit vrel , so dass sich f¨  vrel ρg dp FW (Rep ) = FW (Rep ) f (vrel ) dvrel , Rep = . (6.7) µg R3 63

6 Konstruktion eines adaptiven Simulationsverfahrens vrel,z

FW,z Stokes-Bereich

vrel,z

F W,z FW,x

vrel,x

vrel,x

F W,x FW,z

¨ Ubergangsund Newton -Bereich

FW,z (Rep ) FW,z (Rep ) FW,x FW,x (Rep ) FW,x (Rep )

¨ Abbildung 6.5: Anisotrope Ubertragung der Schwankungen der Relativgeschwindigkeit auf die Widerstandskraft. Aufgrund der anisotropen Schwankungs¨ ubertragung von vrel auf FW sind die Gr¨oßen der Widerstandskraft, wie sie sich durch die Anwendung der beiden Einzelverfahren ergeben, voneinander verschieden. In der Abbildung 6.5 ist diese Problematik f¨ ur einen zweidimensionalen Fall noch einmal schematisch dargestellt. F¨ ur den Stokes-Bereich f¨ uhrt die Berechnung der Widerstandskraft sowohl mit dem Zwei-Kontinua- als auch mit dem Einzelpartikel-Verfahren aufgrund der direkten Proportionalit¨at zwischen σF2 W und σv2 rel zu ¨ und Newton-Bereich ergeben sich den gleichen Ergebnissen (F W,z , F W,x). Im Ubergangszwei unterschiedliche Werte f¨ ur die Widerstandskraft. Im Einzelpartikel-Verfahren f¨ uhrt die Nicht-Linearit¨at zwischen σF2 W und σv2 rel bzw. vrel zu einer Verschiebung der mittleren Widerstandskraft FW (Rep ) zu h¨oheren Werten hin. Der Grund daf¨ ur ist, dass Partikel mit einer h¨oheren Relativgeschwindigkeit vrel einer u ¨berproportional h¨oheren Widerstandskraft ausgesetzt sind, als dies bei Partikel mit geringeren Relativgeschwindigkeiten der Fall ist. Die Berechnung der Widerstandskraft mit dem Zwei-Kontinua-Verfahren ist daher durch die Verwendung von mittleren Werten f¨ ur die Relativgeschwindigkeit mit einem systematischen Fehler behaftet. Um diese Abweichung quantitativ zu erfassen, werden die Widerstandsbeiwerte cd , wie sie sich aus den beiden numerischen Verfahren in Abh¨angigkeit der mittleren Reynolds-Zahl der Partikel Rep und der Standardabweichung der Relativgeschwindigkeit σvrel ergeben, miteinander verglichen. Hierzu wird f¨ ur den Widerstandsbeiwert die in Kapitel 3 eingef¨ uhrte Korrelation nach Turton&Levenspiel verwendet. Der aus dem Einzelpartikel-Verfahren resultierende mittlere Widerstandsbeiwert errechnet sich aus:  cd (Rep ) = cd (Rep ) f (vrel ) dvrel , (6.8) R3

64

6.2 Informationsaustausch an den Zonengrenzen mit der Verteilungsfunktion der Relativgeschwindigkeit:   1 (vrel − vrel )2 f (vrel ) = .

3 exp − 2 σv2 rel 2 π σ2 2

(6.9)

vrel

F¨ ur einen Vergleich zwischen cd (Rep ) und dem mittleren Widerstandsbeiwert, wie er sich aus dem Zwei-Kontinua-Verfahren in Abh¨angigkeit der mittleren Relativgeschwindigkeit vrel ergibt, muss (6.8) numerisch approximiert werden. Bildet man nun aus dem Ergebnis f¨ ur cd (Rep ) und aus dem mittleren Widerstandsbeiwert cd (Rep ) einen Quotienten, so l¨asst sich folgender Korrekturfaktor ckt definieren: ckt =

cd (Rep ) . cd (Rep )

(6.10)

In der Abbildung 6.6 ist ckt in Abh¨angigkeit der mittleren Reynoldszahl der Partikel Rep und der Standardabweichung σvrel dargestellt. F¨ ur erh¨ohte Werte von Rep und σvrel erkennt

1.5 1.4 ckt

1.3 1.2 1.1 1 300

200 100 Rep

0

0

0.6 0.4 0.2 σvrel

0.8

Abbildung 6.6: Korrekturfaktor ckt . man einen erheblichen Anstieg des Korrekturfaktors. In diesem Fall wird die mit dem ZweiKontinua-Verfahren ermittelte Gr¨oße der Widerstandskraft untersch¨atzt. Um diesen offensichtlichen Modellfehler des Zwei-Kontinua-Verfahrens zu beheben, wird der aus der mittleren Relativgeschwindigkeit resultierende, mittlere Widerstandsbeiwert cd (Rep ) mit dem Korrekturfaktor ckt multipliziert. Hierzu soll der in Abbildung 6.6 dargestellte Zusammenhang zwischen ckt und Rep bzw. σvrel durch eine Korrelation wiedergegeben werden. F¨ ur ckt wird folgender Ansatz verwendet:

ckt ≈ 1 + a log b Rep + 1 . (6.11) F¨ ur die Koeffizienten a und b erh¨alt man durch Regression: a = −0.0087 σvrel + 0.2497 σv2 rel − 0.0635 σv3 rel ,

(6.12) 65

6 Konstruktion eines adaptiven Simulationsverfahrens b = 1.6007 − 0.2457 σvrel + 2.5786 σv2 rel − 4.6241 σv3 rel + 2.0599 σv4 rel .

(6.13)

Durch die Ber¨ ucksichtigung der vorstehenden Korrelation ist sichergestellt, dass die mit dem Zwei-Kontinua- sowie mit dem Einzelpartikel-Verfahren ermittelten Werte der Widerstandskraft nahezu identisch sind.

6.2.2

Modellierung der Geschwindigkeitskovarianz

Das Ausmaß der Schwankungsgeschwindigkeit eines Partikelkollektivs wird durch die Angabe der granularen Temperatur Θ beschrieben. W¨ahrend im Einzelpartikel-Verfahren die Berechnung der granularen Temperatur direkt aus den Schwankungsgeschwindigkeiten der Partikel Θ = |vp
|2 /3 ermittelt wird, ist im Falle des Zwei-Kontinua-Verfahrens eine Transportgleichung f¨ ur Θ zu l¨osen (vgl. (3.25)). Wie in Kapitel 3 beschrieben, ber¨ ucksichtigt diese neben dem konvektiven und diffusiven Transport der granularen Temperatur auch die Produktion bzw. Dissipation der Fluktuationsenergie, die aus dem Einfluss der Reibungskraft zwischen den Partikeln und dem umgebenden Gas resultieren. Sie werden in der Gr¨oße Ψ zusammengefasst: Ψp = β(Kgp − 3Θ).

(6.14)

Kgp ist die Geschwindigkeitskovarianz, sie stellt eine Korrelation zwischen den Fluktuationsgeschwindigkeiten vp
einzelner Partikel und den Fluktuationsgeschwindigkeiten der Gasphase u
g dar: Kgp = u
g vp
.

(6.15)

Zur Modellierung von (6.15) wird ein vereinfachtes Kr¨aftegleichgewicht am Einzelpartikel betrachtet [SM98]: 3 ρg dvp = cd (Rep ) |up − vp | (ug − vp ) . dt 4 ρp dp

(6.16)

Nach der Subtraktion der entsprechenden zeitgemittelten Gr¨oßen erh¨alt man die Bewegungsgleichung eines Partikels unter dem Einfluss der Fluktuationsgeschwindigkeit der Gasphase:

dvp
3 ρg = cd (Rep ) |up − vp | u
g − vp
. dt 4 ρp dp

(6.17)

Setzt man nun die G¨ ultigkeit des Stokes-Bereiches f¨ ur den Widerstandsbeiwert cd voraus, so gilt:

dvp
1
u − vp
, = dt τp g

τp =

18 µg , ρp d2p

(6.18)

wobei τp die aerodynamische Partikelrelaxationszeit bezeichnet. Zur L¨osung von (6.18) wird angenommen, dass die von einem Partikel entlang seiner Bahnkurve erfahrene Fluktuation 66

6.2 Informationsaustausch an den Zonengrenzen der Gasphase als ein periodischer Vorgang beschrieben werden kann. So ist es m¨oglich die Fluktuationsgeschwindigkeit der kontinuierlichen Phase durch eine Sinusfunktion mit einer Amplitude A und der Frequenz ω auszudr¨ ucken: u
g (t) = A sin (ω t) .

(6.19)

F¨ ur die eingeschwungene L¨osung der Bewegungsgleichung eines Partikels unter Einwirkung dieser Gasfluktuation gilt: vp
(t) =

u
g (t) ω τp −A cos (ω t) . 1 + ω 2 τp2 1 + ω 2 τp2

(6.20)

Die Bildung des Produktes u
g vp
und eine anschließende zeitliche Mittelung f¨ uhrt u ¨ber einen Koeffizientenvergleich zu folgendem Ergebnis: √ √ u
g vp
= 2k 3Θ. (6.21)

f (Rep )

Dieser Ansatz, der bereits in der Arbeit von Sinclair & Mallo verwendet wurde [SM98], stellt das Produkt der gemittelten Betr¨age der Fluktuationsgeschwindigkeiten dar und darf zun¨achst bei geringen mittleren Reynolds-Zahlen der Partikel Rep