Document [original]

Untersuchung der Dynamik fluider Partikel

auf Basis der Volume of Fluid Methode

Von der Fakultät für Naturwissenschaften

- Department Chemie -

der Universität Paderborn

zur Erlangung des akademischen Grades

Doktor der Naturwissenschaften

- Dr. rer. nat. -

genehmigte

Dissertation

von

Martin Schmidtke

aus Paderborn

Paderborn, Januar 2008

Die vorliegende Arbeit entstand am Institut für Technische Chemie und

Chemische Verfahrenstechnik des Departments Chemie der Universität

Paderborn im Zeitraum von Juni 2003 bis Dezember 2007.

Tag der mündlichen Prüfung: 25.1.2008.

Referent: Professor Dr.-Ing. H.-J.Warnecke

Universität Paderborn

Department Chemie

Institut für Technische Chemie und chemische Verfahrenstechnik

Korreferent: Professor Dr. rer. nat. D. Bothe

RWTH Aachen

Lehrstuhl für Mathematik

Center for Computational Engineering Science

Im Rahmen dieser Arbeit entstanden folgende Veröffentlichungen:

D. Bothe, M. Schmidtke, H.-J. Warnecke:

Direct Numerical Computation of the Lift Force acting on Single Bubbles.

International Conference on Multiphase Flow, ICMF2007, 9.-13. Juli

Leipzig, 2007.

D. Bothe, M. Schmidtke, H.-J. Warnecke:

VoF-simulation of the lift force for single bubbles in a simple shear flow.

Chem. Eng. Tech. 29(9), 1048-1053, 2006.

M. Schmidtke, D. Bothe, H.-J. Warnecke:

VoF-simulation of the rise behaviour of single air bubbles in linear shear flows.

in Proc. of 3

Int. Berlin Workshop on "Transport Phenomena with moving

boundaries" (M. Kraume, Ed.), Berlin 2005.

Danksagung

Viele Menschen haben dazu beigetragen, daß diese Dissertation entstehen

konnte.

Mein besonderer Dank gilt

Professor Dr. Hans-Joachim Warnecke für die Möglichkeit der Mitarbeit in

seinem Arbeitskreis. Durch ihn erhielt ich viele interessante Einblicke in die

technische Chemie. Er hat diese Arbeit in vielfältiger Weise unterstützt, nicht

zuletzt auch in organisatorischer und menschlicher Hinsicht.

Professor Dr. Dieter Bothe, der diese Arbeit mit großem Engagement fachlich

betreut hat. Im interdisziplinären Gespräch mit ihm und Professor Dr. Hans-

Joachim Warnecke entwickelten sich immer wieder neue Blickweisen auf die

Problemstellungen, die zu weiteren Untersuchungen anspornten.

Professor Dr. Marsmann und Professor Dr. Huber für ihre Bereitschaft, als

Mitglieder der Prüfungskommission zu wirken.

Dr. Kerstin Wielage und Michael Kröger, die den Code FS3D gewartet und für

diese Arbeit erweitert haben, sowie Dr. Mario Koebe, der mich in die

Benutzung von FS3D eingewiesen hat.

Dem Paderborn Center for Parallel Computing (PC

) für die zur Verfügung

gestellte Rechenleistung und die gute Benutzerberatung seitens der

Systemadministratoren, insbesondere durch Axel Keller.

Meinen ehemaligen Kollegen Bhaskar Bandarapu und Houman Shirzadi im

„Multikulti-Büro“, sowie allen Mitarbeitern der technischen Chemie für die gute

und fruchtbare Arbeitsatmosphäre.

Meinen neuen Kollegen am Forschungszentrum Dresden-Rossendorf für ihre

technische und organisatorische Unterstützung beim Aufschreiben dieser

Dissertation.

Dr. Alexander Grahn für nützliche Hinweise im Umgang mit LATEX.

Schließlich möchte ich in besonderer Weise meinen Eltern danken, die mich

während dieser Arbeit nach Kräften unterstützt haben. Als fleißiger Lektor

‚promovierte’ mein Vater verdientermaßen zum „Doktor-Vater“.

Inhaltsverzeichnis

1. Einleitung 1

2. Hydrodynamik fester und fluider Partikel 5

2.1. Der Str¨

omungswiderstand fester Kugeln . . . . . . . . . . . . . . . . 5

2.2. Das Diagramm von Clift, Grace und Weber . . . . . . . . . . . . . . 6

3. Kontinuumsmechanische Modellierung 9

3.1. Bilanzierung von Zweiphasenstr¨

omungen . . . . . . . . . . . . . . . . 9

3.1.1. Integrale Bilanzen . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.1.2. Differenzielle Bilanzen . . . . . . . . . . . . . . . . . . . . . . 13

3.1.3. Sprungbedingungen . . . . . . . . . . . . . . . . . . . . . . . . 14

3.1.4. Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . 15

4. Numerische Verfahren 17

4.1. Die VoF-Methode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

4.2. Implementierung der Grenzfl¨

achenkraft . . . . . . . . . . . . . . . . . 19

5. Dimensionsanalyse der Navier-Stokes-Gleichungen 21

5.1. Die einphasigen Navier-Stokes-Gleichungen in dimensionsloser Form . 21

5.2. Die mehrphasigen Navier-Stokes-Gleichungen in dimensionsloser Form 24

6. Validierung des Codes FS3D 27

6.1. Die Hadamard-Rybczinski-L¨

osung . . . . . . . . . . . . . . . . . . . . 27

6.2. Numerische L¨

osungen f¨

ur langsame, sph¨

arische fluide Partikel . . . . 29

I. Bewegung von ¨

Oltropfen in Wasser 37

Inhaltsverzeichnis

7. Simulation von Mais¨

oltropfen 41

7.1. Konfiguration des Rechengebietes . . . . . . . . . . . . . . . . . . . . 41

7.2. Einfluß des Wandabstandes auf die Aufstiegsgeschwindigkeit . . . . . 41

7.3. Geschwindigkeiten innerhalb des Tropfens . . . . . . . . . . . . . . . 42

7.4. Trajektorien . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

7.5. Aufstiegsgeschwindigkeiten und Tropfenformen . . . . . . . . . . . . . 51

7.6. Einfluß der tr¨

agen Masse fluider Partikel auf das Aufstiegsverhalten . 52

8. Rizinus¨

ol 55

8.1. Aufstiegsgeschwindigkeiten und Tropfenformen . . . . . . . . . . . . . 55

9. Widerstandsbeiwerte f¨

ur ¨

Oltropfen 57

10.Zusammenfassung des ersten Teils 59

II. Blasenaufstieg in linearen Scherstr¨

omungen 61

11.Stand des Wissens 63

11.1. Blasen in vertikalen Rohren . . . . . . . . . . . . . . . . . . . . . . . 63

11.2. Theoretische Untersuchungen zu kugelf¨

ormigen Blasen in linearen Scher-

str¨

omungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

11.3. Das Experiment von Tomiyama . . . . . . . . . . . . . . . . . . . . . 67

12.Simulation von Blasen in Glycerinl¨

osungen 69

12.1. Randbedingungen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

12.2. Blasenformen und Aufstiegsgeschwindigkeiten . . . . . . . . . . . . . 70

12.3. Aufstiegsbahnen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

12.4. Die Bestimmung des Liftkoeffzienten CLaus Blasentrajektorien . . . 76

12.5. Die Abh¨

angigkeit des Liftkoeffizienten von der Scherrate . . . . . . . 78

12.6. CLin Abh¨

angigkeit vom Blasendurchmesser . . . . . . . . . . . . . . 79

12.7. Abh¨

angigkeit der simulierten Liftkoeffizienten von der Aufl¨

osung . . . 80

12.8. Variation der Oberfl¨

achenspannung . . . . . . . . . . . . . . . . . . . 83

12.9. Variation der Gravitation . . . . . . . . . . . . . . . . . . . . . . . . . 87

12.10. Luftblasen in Wasser . . . . . . . . . . . . . . . . . . . . . . . . . . 88

13.Untersuchung der Geschwindigkeits- und Druckfelder an der Phasen-

grenze 93

13.1. Zweidimensionale Blasen . . . . . . . . . . . . . . . . . . . . . . . . . 93

13.1.1. Zirkulation um zweidimensionale Blasen . . . . . . . . . . . . 94

13.1.2. Dynamischer Druck um zweidimensionale Blasen . . . . . . . 95

13.2. Dreidimensionale Blasen . . . . . . . . . . . . . . . . . . . . . . . . . 98

13.2.1. Zirkulation um dreidimensionale Blasen . . . . . . . . . . . . . 98

13.2.2. Dynamischer Druck und Scherkraft an der Phasengrenze . . . 101

13.2.3. Validierung des Verfahrens zur Berechnung der Kr¨

afte an der

Phasengrenze . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

14.Horizontalgeschwindigkeiten im Blasennachlauf 107

15.Zusammenfassung des zweiten Teils 113

A. Anhang 115

A.1. Symbolverzeichnis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

A.2. Das Druckfeld in einem Hill’schen Wirbel . . . . . . . . . . . . . . . . 117

A.3. Konstruktion einer einh¨

ullenden Fl¨

ache parallel zur Phasengrenze . . 120

A.4. Geschwindigkeitsfeld um einen oszillierenden ¨

Oltropfen . . . . . . . . 122

Literaturverzeichnis 125

Inhaltsverzeichnis

1. Einleitung

In vielen nat¨

urlichen oder verfahrenstechnischen Prozessen treten Str¨

omungen auf, in

denen zwei Fluide gemeinsam transportiert werden. Ist die gegenseitige L¨

oslichkeit

gering, so bildet sich zwischen den beiden Fluiden eine Grenzschicht, die sogenann-

te Phasengrenze. In diesem Falle spricht man von einer zweiphasigen Str¨

omung. Die

Fluide werden auch Phasen genannt. Zweiphasige Str¨

omungen mit einer beweglichen

Phasengrenze bestehen entweder aus zwei fl¨

ussigen Phasen (etwa Wasser und ¨

Ol)

oder einer fl¨

ussigen und einer gasf¨

ormigen Phase (z.B. Wasser und Luft). Zwei Gase

hingegen k¨

onnen sich auf molekularer Ebene vermischen. Hier tritt keine lokalisierba-

re Phasengrenze auf, weshalb sie einphasige Systeme bilden. Ein anderer Typus von

Zweiphasenstr¨

omungen ist der Transport fester Partikel in einem Fluid, beispielswei-

se von Sand in Wasser oder in Luft. Dabei wird die Gesamtheit der Partikel als eine

(feste) Phase bezeichnet. Die Oberfl¨

ache der einzelnen Partikel stellt in diesem Fall

eine starre Phasengrenze dar.

Wenn eine Phase auf getrennte Gebiete verteilt ist, die von der zweiten zusam-

menh¨

angenden Phase vollst¨

andig umschlossen sind, so wird die zusammenh¨

angende

Phase kontinuierliche Phase genannt. Die andere Phase wird als disperse Phase be-

zeichnet (siehe Abbildung 1.1 auf Seite 3). Die Gebiete mit disperser Phase bilden

feste oder fluide Partikel. Fluide Partikel sind je nach Aggregatzustand der dispersen

Phase entweder Tropfen oder Blasen. In der vorliegenden Arbeit wird das Bewe-

gungsverhalten einzelner fluider Partikel mittels numerischer Simulation untersucht.

In Teil 1 wird der Aufstieg von ¨

Oltropfen in Wasser behandelt. In Teil 2 geht es um

die Bewegung von Luftblasen in Wasser oder w¨

assrigen L¨

osungen. Die Motivation f¨

diese Untersuchungen und die entsprechenden verfahrenstechnischen Anwendungen

sind zu Beginn der jeweiligen Teile skizziert.

Mittels der in dieser Arbeit verwendeten Volume-of-Fluid (VoF) Methode ist es

m¨

oglich, die Bewegung der Phasengrenze und damit die Form¨

anderungen eines flui-

den Partikels zu erfassen. Der numerische Aufwand ist recht hoch, weshalb bei der

1. Einleitung

heute zur Verf¨

ugung stehenden Rechnerleistung mit dieser Methode eine gemeinsame

Bewegung f¨

ur h¨

ochstens einige wenige Partikel simuliert werden kann. In technisch re-

levanten Prozessen hingegen bewegen sich tausende Tropfen oder Blasen gleichzeitig.

Will man derartige Prozesse unter Verwendung von Simulationen optimieren, werden

daher vereinfachende Modelle verwendet, welche die individuelle ver¨

anderliche Form

der einzelnen Partikel vernachl¨

assigen. So werden beispielsweise in der Euler-Lagrage-

Methode die Partikel als punktf¨

ormige Objekte betrachtet, an denen verschiedene

Kr¨

afte als sogenannte Punktkr¨

afte ansetzen. Die Bahn der einzelnen Partikel wird

berechnet, ebenso wie die Entwicklung der Str¨

omung der kontinuierlichen Phase. Ei-

ne weitere Vergr¨

oberung der Betrachtungsweise f¨

uhrt zur Euler-Euler-Methode. Hier

wird die Gesamtheit der festen oder fl¨

ussigen Partikel als Fluid betrachtet, welche

das Fluid der kontinuierlichen Phase vollst¨

andig durchdringt. Die Individualit¨

at der

Partikel wird dabei vern¨

achl¨

assigt. Die Verteilung der Partikel wird nur noch in Form

einer ¨

ortlichen Partikeldichte beschrieben, beispielsweise als Skalarfeld αp, welches

f¨

ur jeden Ort den Volumenanteil der dispersen Phase angibt. Zur Schließung abgelei-

teter Modelle ist die Beschreibung der Kr¨

afte auf disperse Partikel in Abh¨

angigkeit

von lokalen Str¨

omungseigenschaften sowie von Gr¨

oße und Form der Partikel zwingend

erforderlich.

Derartige Partikelkr¨

afte sind beispielsweise die Auftriebskraft und die Widerstands-

kraft, welche der Bewegung des Partikels durch die kontinuierliche Phase entgegen-

wirkt. Der Impulsaustausch zwischen kontinuierlicher und disperser Phase l¨

aßt durch

die Partikelkr¨

afte beschreiben. Wird beispielsweise einerseits die Bewegung der Par-

tikel durch die Widerstandskraft abgebremst, so wird andererseits (nach dem Prinzip

”actio gleich reactio“) gleichzeitig das die Partikel umgebende kontinuierliche Fluid

beschleunigt. Die Genauigkeit der f¨

ur die Auslegung technischer Apparaturen benutz-

ten Simulationsmethoden h¨

angt entscheidend davon ab, wie akkurat die verwendeten

Modelle f¨

ur die an den fluiden Partikeln wirksamen Kr¨

afte sind.

Die in dieser Arbeit vorgestellten Untersuchungen dienen einerseits dazu, mehr Ver-

st¨

andnis ¨

uber die Entstehung der an fluiden Partikeln wirkenden Kr¨

afte zu gewin-

nen; andererseits dazu, Modelle f¨

ur Punktkr¨

afte zu verbessern oder zu best¨

atigen. Im

Teil 1 geht es in erster Linie um die Widerstandskraft an in Wasser aufsteigenden

Oltropfen und der sich daraus ergebenden Aufstiegsgeschwindigkeit. Teil 2 behan-

delt neben der Widerstandskraft an Blasen haupts¨

achlich die sogenannte Liftkraft,

die in Scherstr¨

omungen auftritt. Diesen beiden Hauptteilen sind einige einf¨

uhrende

Kapitel vorangestellt, in denen die relevanten theoretischen Grundlagen, numerischen

Methoden und empirischen Korrelationen beschrieben sind.

00000000

11111111

000

111

0000000

1111111

0000

1111

0000

1111

Abbildung 1.1.: Verteilung von kontinuierlicher und disperser Phase in einem Gebiet.

Blau: Kontinuierliche Phase. Weiß (schraffiert): Disperse Phase.

1. Einleitung

2. Hydrodynamik fester und fluider Partikel

2.1. Der Str¨

omungswiderstand fester Kugeln

Weil die Kugel der einfachste geometrische K¨

orper ist, beziehen sich die meisten

Untersuchungen zu Partikelbewegungen in Fl¨

ussigkeiten auf Kugeln. Der Str¨

omungs-

widerstand eines K¨

orpers kann durch eine dimensionslose Kennzahl CDcharakteri-

siert werden. Diese Kennzahl CDist im Allgemeinen formspezifisch und wird Wider-

standsbeiwert (engl.: drag coefficient) genannt. Sie ist definiert als Quotient aus der

Widerstandskraft FWund der Kraft des dynamischen Druckes auf der angestr¨

omten

Querschnittsfl¨

ache Ades K¨

orpers:

CD=FW

2ρw2A.(2.1)

Hier ist ρdie Dichte des Fluids und wdie Str¨

omungsgeschwindigkeit. Die angestr¨

omte

Querschnittsfl¨

ache Aist die Fl¨

ache der Projektion des K¨

orpers in Str¨

omungsrichtung.

F¨

ur eine Kugel ist dies eine Kreisfl¨

ache. Weil Aim Nenner von (2.1) auftritt, ist der

Widerstandsbeiwert nicht von der Gr¨

oße eines K¨

orpers abh¨

angig, sondern ausschließ-

lich von der Form. Es besteht jedoch eine Abh¨

angigkeit des Widerstandsbeiwertes von

der Geschwindigkeit w, beziehungsweise von der zu wproportionalen Reynoldszahl.

Die Reynoldszahl ist eine dimensionslose Kennzahl, die das Verh¨

altnis von Tr¨

agheits-

kraft zu Z¨

ahigkeitskraft angibt,

Re = ρLw

µ.(2.2)

Hier ist µdie Viskosit¨

at des Fluids und Leine charakteristische L¨

ange des K¨

orpers.

Im Falle einer Kugel ist der Durchmesser die charakteristische L¨

ange. F¨

ur sehr klei-

ne Str¨

omungsgeschwindigkeiten (Re≪1) kann der Widerstandsbeiwert von Kugeln

analytisch hergeleitet werden (vgl. Abschnitt 6.1).

CD=24

Re f¨

ur Re ≪1.(2.3)

2. Hydrodynamik fester und fluider Partikel

F¨

ur gr¨

oßere Reynoldszahlen kann der Widerstandsbeiwert nur experimentell ermit-

telt werden. Eine besonders einfache Methode ist es, eine Kugel so lange im fluiden

Medium aufsteigen oder fallen zu lassen, bis sie eine konstante Endgeschwindigkeit

wTerreicht hat. Dann sind Auftriebskraft und Widerstandskraft im Gleichgewicht:

2ρ wT2A=g∆ρV . (2.4)

Hier ist ∆ρdie Differenz zwischen Kugeldichte und der Dichte des fluiden Mediums.

Setzt man f¨

ur die Querschnittsfl¨

ache A=πd2/4 und f¨

ur das Volumen V=πd3/6 der

Kugel in (2.4) ein, so ergibt sich mit dals Kugeldurchmesser

CD=4

∆ρ

.(2.5)

Es gibt zahlreiche Korrelationen, die einen unterschiedlichen Geltungsbereich in Be-

zug auf die Reynoldszahlen haben. Einen sehr großen Geltungsbereich deckt die Kor-

relation von Turton und Levenspiel ab, die auf einer sehr breiten Datenbasis vieler

verschiedener Experimentatoren beruht (Turton und Levenspiel, 1986). Danach ist

CD=24

Re 1 + 0,173 Re0,675+0,413

1 + 16300Re−1,09 f¨

ur Re <105.(2.6)

Im Intervall 103<Re <105ist der Widerstandsbeiwert nahezu konstant (CD≈0,47 ,

siehe Abbildung 7.2, S. 44). Dies wurde bereits von Newton beobachtet, weshalb man

hier auch von dem Newton’schen Bereich spricht. Newton hatte hierzu Fallexperi-

mente von der Kuppel der St.-Pauls-Kathedrale in London durchgef¨

uhrt.

Wenn feste Kugeln in Fl¨

ussigkeiten frei aufsteigen oder sinken, kann durch eine peri-

odische Wirbelabl¨

osung im Nachlauf eine zickzack-f¨

ormige Bewegung auftreten (Hal-

jasmaa, 2006). Mit dem Einsetzen einer derartigen seitlichen Oszillation ist eine deut-

liche Absenkung der Vertikalgeschwindigkeit verbunden. Demzufolge liegt der Wider-

standsbeiwert f¨

ur oszillierende Kugeln deutlich h¨

oher als nach (2.6). Eine seitliche

Oszillation der Kugel kann unterbunden werden, indem man sie in einem Str¨

omungs-

kanal fixiert. Gleichung (2.6) gilt f¨

ur fixierte Kugeln sowie f¨

ur geradlinig aufsteigende

oder sinkende Kugeln.

2.2. Das Diagramm von Clift, Grace und Weber

Im Gegensatz zu festen Kugeln besitzen fluide Partikel eine bewegliche Phasengren-

ze. Dadurch sind Formver¨

anderungen m¨

oglich. Die Partikelform kann grob durch das

2.2. Das Diagramm von Clift, Grace und Weber

Durchmesserverh¨

altnis dV/dHcharakterisiert werden, wobei dVder vertikale und dH

der horizontale Durchmesser sind. Die Partikelgr¨

oße wird entweder durch sein Volu-

men Vbeschrieben oder durch den ¨

Aquivalenzdurchmesser d, der den Durchmesser

einer Kugel gleichen Inhalts angibt:

d= (6V/π)1/3.(2.7)

Je nach Gr¨

oße und Aufstiegsgeschwindigkeit werden f¨

ur fluide Partikel vier verschiede-

ne Formklassen unterschieden: Kugelf¨

ormig (sehr klein oder sehr langsam), ellipsoid,

wobbelnd (d.h. amorph, formver¨

anderlich) und kappenf¨

ormig. Die Kappenform wird

nur bei sehr großen Blasen beobachtet. Die gr¨

oßten Tropfen sind wobbelnd. Tropfen

in der Gr¨

oße von Kappenblasen sind nicht stabil; sie zerfallen. F¨

ur kugelf¨

ormige flui-

de Partikel beobachtet man in Experimenten oft die gleiche Aufstiegsgeschwindigkeit

bzw. Widerstandsbeiwerte wie bei festen Kugeln. Offenbar findet eine Bewegung der

inneren Grenzschicht parallel zur Phasengrenze nicht statt. Diese Immobilisierung

der Phasengrenze kann durch Anlagerung von oberfl¨

achenaktiven Molek¨

ulen (soge-

nannten ”Surfactants“) verursacht werden. Derartige Stoffe (wie z.B. Tenside) sind

in Fl¨

ussigkeiten stets enthalten, wenn sie nicht aufwendig gereinigt werden.

Die Formenvielfalt fluider Partikel macht es außerordentlich schwierig, universelle

Korrelationen f¨

ur ihr Aufstiegsverhalten zu finden. Daher sind Korrelationen f¨

ur die

Aufstiegsgeschwindigkeit oder die Partikelform im Allgemeinen auf einen engen Wer-

tebereich bez¨

uglich der Partikelgr¨

oße oder der Stoffeigenschaften beschr¨

ankt. Eine

Ausnahme bildet hier das Diagramm von Clift, Grace und Weber, welches f¨

ur einen

außerordentlich großen Wertebereich eine M¨

oglichkeit zur Absch¨

atzung der Formklas-

se und der Aufstiegsgeschwindigkeit fluider Partikel bietet. Die Aufstiegsgeschwindig-

keit wird in Form der Reynoldszahl Rec=ρcwTd/µcdes Partikels (bezogen auf das

umgebende kontinuierliche Fluid mit der Dichte ρcund der Viskosit¨

at µc) angege-

ben. Das Diagramm 2.1 gibt Recin Abh¨

angigkeit von zwei weiteren dimensionslosen

Kennzahlen an, die ihrerseits von der Partikelgr¨

oße und den Stoffeigenschaften der

Fluide abh¨

angen. Diese sind die E¨

otv¨

oszahl

E¨

o = g∆ρ d2

σ,(2.8)

welche die Auftriebskraft zur Oberfl¨

achenspannung ins Verh¨

altnis setzt, und die

2. Hydrodynamik fester und fluider Partikel

Morton-Zahl

Mo = g∆ρµ4

ρ2

cσ3,(2.9)

die die Stoffeigenschaften charakterisiert. Ausgenommen vom Geltungsbereich dieser

Korrelation von Clift et al. (1978) sind Tropfen in Gas (etwa Regentropfen) d.h. Situa-

tionen, in denen Dichte und Viskosit¨

at des Partikels um mehr als eine Gr¨

oßenordung

h¨

oher ist als die des kontinuierlichen Fluids.

Abbildung 2.1.: Diagramm zur Absch¨

atzung von Form und Geschwindigkeit fluider Par-

tikel von Clift et al. (1978).

3. Kontinuumsmechanische Modellierung

3.1. Bilanzierung von Zweiphasenstr¨

omungen

Die hier verwendete kontinuumsmechanische Modellierung setzt voraus, das die be-

trachteten Volumenelemente eine große Anzahl von Molek¨

ulen enthalten, damit in-

nerhalb jeder Phase die makroskopischen Gr¨

oßen (wie zum Beispiel die Dichte ρ)

als stetige Funktionen beschrieben werden k¨

onnen. Die Phasengrenze Γ zeichnet sich

dadurch aus, daß mindestens eine dieser Gr¨

oßen an Γ unstetig ist. Im Allgemeinen

unterscheiden sich zwei Fluide sowohl in der Dichte als auch in der Viskosit¨

at, so-

daß beide Gr¨

oßen an der Phasengrenze unstetig sind. Die Phasengrenze selbst wird

als masselose st¨

uckweise glatte Fl¨

ache angenommen. Im Folgenden werden zun¨

achst

die integralen und dann die differenziellen Bilanzen f¨

ur Zweiphasenstr¨

omungen auf-

gestellt.

3.1.1. Integrale Bilanzen

Aus der Bilanzierung der Gr¨

oßen Masse, Impuls und Energie ¨

uber ein Kontrollvolu-

men, welches die Phasengrenze schneidet, ergeben sich die integralen Bilanzgleichun-

gen f¨

ur zweiphasige Str¨

omungen.

Zwei auf molekularer Ebene nicht vermischbare (d.h. gegenseitig nicht l¨

osbare) Fluide

befinden sich jeweils in den Teilgebieten Ω−(t) und Ω+(t) im Gebiet G⊂R3. Die

Phasengrenze Γ(t) zwischen den Teilgebieten bewegt sich mit der Geschwindigkeit uΓ.

Die Fl¨

achennormale nΓzeigt in Richtung Ω+(t). Das f¨

ur die Bilanzierung verwendete

Kontrollvolumen schneidet die Phasengrenze.

3. Kontinuumsmechanische Modellierung

Massenbilanz

Ein ortsfestes Volumen Venth¨

alt die Masse RVρdV . Wird ein Fl¨

achenelement dA

auf seinem Rand mit der Geschwindigkeit udurchstr¨

omt, so fließt durch dA pro

Zeiteinheit die Masse ρu·ndA. Aus der Massenerhaltung folgt: Die Massen¨

anderung

in Vist gleich der Masse, die ¨

uber den Rand ∂V fließt. Also gilt

dt ZV

ρ dV =−Z∂V

ρu·ndA , (3.1)

oder umgestellt

dt ZV

ρ dV +Z∂V

J·ndA = 0 mit J=ρu.(3.2)

Die vektorielle Gr¨

oße Jwird Massenstromdichte genannt. F¨

ur eine Fl¨

ache mit der

Normalen ngibt J·ndie Rate des Massenflusses pro Fl¨

acheneinheit an. Im Sonder-

fall, daß ρinnerhalb jeder Phase konstant ist, tritt eine Massen¨

anderung im Kon-

trollvolumen Vgenau dann ein, wenn sich die Phasengrenze durch Vbewegt. Die

Dichte ρan einem Punkt innerhalb von Vh¨

angt dann davon ab, auf welcher Seite

der Phasengrenze sich der Punkt befindet.

Impulsbilanz

Es wird hier ein Volumen V(t) betrachtet, welches mit der Str¨

omung mitbewegt wird.

Der Impuls dieses Volumens betr¨

agt

I(t) = ZV(t)

ρudV . (3.3)

Die Rate der Impuls¨

anderung eines K¨

orpers ist nach dem 2. Newton’schen Gesetz die

Summe aller angreifenden Kr¨

afte

dtI(t) = X

Fi.(3.4)

Es wird grunds¨

atzlich zwischen zwei verschiedenen Krafttypen unterschieden:

K¨

orperkr¨

afte wirken auf alle Teilchen eines Volumens durch ein Kraftfeld. Es kann

sich beispielsweise um ein elektrisches, magnetisches oder Gravitationsfeld handeln.

Die massenspezifische Kraftdichte ρfwirkt im gesamten Volumen. Daher werden sie

auch Volumenkr¨

afte genannt.

3.1. Bilanzierung von Zweiphasenstr¨

omungen

Kontaktkr¨

afte wirken an der Oberfl¨

ache eines Volumens, weshalb sie auch Ober-

fl¨

achenkr¨

afte genannt werden. Ihre fl¨

achenspezifische Dichte S∂V (t)hat die Einheit

Kraft pro Fl¨

ache. Beispiele sind Druckkraft und Schubspannungen.

F¨

ur die Impuls¨

anderung gilt demnach:

dt ZV(t)

ρudV =S∂V (t)+ZV(t)

ρfdV . (3.5)

In idealen (d.h. reibungsfreien) einphasigen Fluiden treten keine Schubspannungen

auf. Als Oberfl¨

achenkraft tritt dann nur der Druck pin Erscheinung. Auf ein Fl¨

achen-

element dA mit der Normalen nwirkt die Kraft pndA, d.h.

S∂V (t)=−Z∂V (t)

pndA . (3.6)

In reibungsbehafteten Fluiden treten Spannungskr¨

afte auf. In zweiphasigen Systemen

kommen Grenzfl¨

achenkr¨

afte hinzu. Dadurch erweitert sich die Impulsbilanz zu

dt ZV(t)

ρudV =ZV(t)

ρfdV +Z∂V (t)

T·ndA +FA(t).(3.7)

Hierbei ist der Spannungstensor Teine 3 ×3-Matrix mit T=−pI+S. Darin ist

Idie Einheitsmatrix und Sder viskose Spannungstensor, welcher die durch innere

Reibung verursachten Spannungskr¨

afte beschreibt. Schneidet das Kontrollvolumen

V(t) die Grenzfl¨

ache Γ(t), so liefert die Schnittfl¨

ache A= Γ(t)∩V(t) den Beitrag

FA(t)durch die Grenzfl¨

achenspannung. Um diesen Beitrag herzuleiten, betrachten wir

den Rand von A, welcher die geschlossene Kurve Cbildet (siehe Abbildung 3.1). Wenn

V(t)

(t)

ΓA

Abbildung 3.1.: Schnitt eines Kontrollvolumens V(t) mit Phasengrenze Γ(t).

3. Kontinuumsmechanische Modellierung

dsein Linienelement von Cund die ¨

außere Fl¨

achennormale an dieser Stelle nΓist, so

wirkt durch die Oberfl¨

achenspannung σauf dieses Element die Kraft dFA=σs×nΓ

tangential zu Γ. Durch eine Integration ¨

uber die Randkurve Cergibt sich als gesamte

Kraft FA(t), welche auf die Schnittfl¨

ache A(t) wirkt:

FA(t) = IC

σ(ds×nΓ) = IC

ds×σnΓ.(3.8)

Unter Verwendung des Stokes’schen Satzes kann das Linienintegral um den Fl¨

achen-

rand Cin ein Fl¨

achenintegral ¨

uber Aumgewandelt werden. Dies liefert

FA=ZA

(nΓ×∇)×σnΓdA . (3.9)

Der Integrand beschreibt die Fl¨

achendichte der Grenzfl¨

achenspannung fΓ. Sie l¨

aßt

sich umformen zu:

fΓ=∇σ−nΓ(nΓ·∇σ) + σnΓ·∇nΓ−σnΓ(∇·nΓ).(3.10)

Der zweite Term stellt den Normalanteil von ∇σdar, so daß die beiden ersten Terme

zusammen den Tangentialanteil von ∇σergeben. Dieser ist der Fl¨

achengradient ∇Γσ.

Der dritte Term ist gleich Null, und so ergibt sich:

fΓ=∇Γσ−σnΓ(∇·nΓ).(3.11)

Wenn die Oberfl¨

achenspannung σkonstant ist, entf¨

allt auch der Term ∇Γσ. Im Allge-

meinen ist σtemperaturabh¨

angig und abh¨

angig von der Konzentration oberfl¨

achen-

aktiver Stoffe (Surfactants). Diese Stoffe senken zumeist die Oberfl¨

achenspannung.

Treten r¨

aumliche Gradienten in der Temperatur oder in der Konzentration der Sur-

factants auf, so ist ∇Γσungleich Null. Die hierdurch bewirkten Ph¨

anomene werden als

Marangoni-Effekte bezeichnet. Der Term −∇·nΓist eine skalare Gr¨

oße und beschreibt

die Kr¨

ummung κΓder Fl¨

ache, die sich aus der Summe ihrer Hauptkr¨

ummungen ergibt

(vgl. Gleichung 3.23). Mit

κΓ=−∇·nΓ(3.12)

erh¨

alt man f¨

ur die Kraft der Oberfl¨

achenspannung auf die Fl¨

ache A:

FΓ=ZA

fΓdA mit fΓ=∇Γσ+σκΓnΓ.(3.13)

F¨

ur ein ortsfestes Volumenelement Vlautet somit die Impulsbilanz:

dt ZV

ρudV +Z∂V

J·ndA =ZV

ρfdV +ZV∩Γ(t)

fΓdA (3.14)

mit dem Fluß J=ρu⊗u−Tund nΓaus (3.13).

3.1. Bilanzierung von Zweiphasenstr¨

omungen

3.1.2. Differenzielle Bilanzen

Die integralen Bilanzen beziehen sich auf ausgedehnte Volumenelemente. Durch den

Grenz¨

ubergang |V| → 0 k¨

onnen sie in differentielle Bilanzen ¨

uberf¨

uhrt werden, die

punktweise gelten. Wird als Bilanzvolumen beispielsweise die Kugel V=Br(x0)

gew¨

ahlt, so geschieht der Grenz¨

ubergang |V| → 0 durch Verkleinern des Kugelradius r

als Grenzwertbildung r→0. Es ergibt sich ein Grenzwert f¨

ur den Quotienten aus dem

Integral der bilanzierten Gr¨

oße ¨

uber das Volumen Vund dem Volumeninhalt |V|.

Massenbilanz

Das Oberfl¨

achenintegral auf der rechten Seite der Gleichung (3.1) kann mittels des

Divergenztheorems in ein Volumenintegral umgewandelt werden. Es ergibt sich die

differentielle Form der Massenbilanz, die auch Kontinuit¨

atsgleichung genannt wird:

∂

∂t ρ+ div (ρu) = 0 .(3.15)

Ist die Dichte ρkonstant, so vereinfacht sich die Kontinuit¨

atsgleichung zu

div u= 0 .(3.16)

Str¨

omungen, f¨

ur die dies gilt, nennt man inkompressibel. In diesem Fall ist (3.15)

aquivalent zu

div u= 0 und ∂

∂t ρ+u·∇ρ= 0 .(3.17)

Dies reduziert sich f¨

ur Fluide konstanter Dichte auf die Gleichung (3.16).

Impulsbilanz

Die differentielle Form der Impulsbilanz ergibt sich, wenn das Divergenztheorem auf

die Oberfl¨

achenintegrale in der integralen Impulsbilanz (3.14) angewandt wird. Wenn

die Gravitation (in Richtung der negativen z-Achse) die einzige Volumenkraft ist, so

gilt

∂

∂t(ρu) + div (ρu⊗u−T) = −ρgez,(3.18)

mit dem Spannungstensor T=−pI+S. F¨

ur Newton’sche Fluide l¨

aßt sich der Z¨

ahig-

keitsspannungstensor Sdurch folgende konstituierende Gleichung ausdr¨

ucken:

S=η(div u)I+µ(∇u+ (∇u)T),(3.19)

3. Kontinuumsmechanische Modellierung

mit der Bulkviskosit¨

at ηund der dynamischen Viskosit¨

at µ. Im Falle eines dichte-

best¨

andigen Fluids vereinfacht sich die Impulsbilanz zu:

∂

∂tu+∇·(u⊗u) = −1

ρ∇p+µ

ρ∆u.(3.20)

Der Quotient µ/ρ =νwird als kinematische Viskosit¨

at bezeichnet. Die Gleichungen

(3.20) und (3.15) bilden zusammen Navier-Stokes-Gleichungen. F¨

ur ideale Fluide,

d.h. solche ohne innere Reibung und damit verbundener Energiedissipation, ergeben

sich aus den Navier-Stokes-Gleichungen die Euler-Gleichungen

∂

∂tu+∇·(u⊗u) = −1

ρ∇pund

div u= 0 .

3.1.3. Sprungbedingungen

Bei Zweiphasenstr¨

omungen sind im Allgemeinen sowohl die Dichte als auch die Vis-

kosit¨

at an der Phasengrenze unstetig. Der Sprung einer Gr¨

oße Ψ an der Phasengrenze

Γ(t) wird als [Ψ] := Ψ+−Ψ−beschrieben, wobei Ψ+und Ψ−die Grenzwerte beider-

seits der Phasengrenze bezeichnen. Wenn kein Phasen¨

ubergang (Verdampfung oder

Kondensation) stattfindet, sind die Normalkomponenten der Geschwindigkeiten an

der Grenzfl¨

ache identisch.

[un] = 0 auf Γ(t).

F¨

ur viskose Systeme wird Haftung an der Phasengrenze angenommen, d.h. die Tan-

gentialkomponente uτist an der Phasengrenze stetig:

[uτ] = 0 auf Γ(t).

Aus der Stetigkeit beider Geschwindigkeitskomponenten auf der Phasengrenze folgt

die Stetigkeit der gesamten Geschwindigkeit an der Phasengrenze. F¨

ur die Massenbi-

lanz lautet die Sprungbedingung demnach

[u] = 0 auf Γ(t).

F¨

ur den Impuls gilt an der Phasengrenze

[ρu⊗(u−uΓ)−T]·nΓ=∇Γσ+σκΓnΓ.(3.21)

3.1. Bilanzierung von Zweiphasenstr¨

omungen

F¨

ur eine r¨

aumlich konstante Grenzfl¨

achenspannung σvereinfacht sich (3.21) zu:

(p+−p−)nΓ−[S]·nΓ=σκΓnΓ.(3.22)

In ruhenden Fluiden ist u= 0 und somit auch S=0. Dann geht Gleichung (3.22) in

die Young-Laplace-Gleichung ¨

uber

(p+−p−) = σκΓ=σ1

r2,(3.23)

mit den Kr¨

ummungsradien r1und r2. Eine ausf¨

uhrliche Herleitung der Sprungbilan-

zen findet sich bei (Ishii, 1975; Slattery, 1999).

3.1.4. Randbedingungen

Um ein System von Differentialgleichungen eindeutig l¨

osen zu k¨

onnen, ist die Kenntnis

der Randbedingungen erforderlich. Wenn ein fester Wert f¨

ur eine Gr¨

oße auf dem Rand

vorgegeben wird, spricht man von einer Dirichlet-Bedingung. Ist die Ableitung am

Rand in normaler Richtung festgelegt, so liegt eine Neumann-Bedingung vor. In der

Str¨

omungssimulation m¨

ussen f¨

ur den Rand des Rechengebietes physikalisch m¨

oglichst

sinnvolle Randbedingungen f¨

ur die Geschwindigkeit eingesetzt werden. Im Folgenden

werden die h¨

aufigsten Randbedingungen aufgef¨

uhrt:

Haftbedingung: Das Str¨

omungsgebiet ist durch eine feste Wand begrenzt. Sie ist

f¨

ur Fluide undurchdringlich, d.h. die Str¨

omungsgeschwindigkeit hat an der Wand

keine normale (zur Wand senkrechte) Komponente un. Durch Reibung haftet eine

Grenzschicht an der Wand. Tangentiale (zur Wand parallele) Geschwindigkeiten uτ

treten demnach ebenfalls nicht auf:

un= 0 und uτ= 0 .

Schlupf-Randbedingung: Der Rand ist f¨

ur das Fluid undurchdringlich. An der Wand

gibt es keine Haftung und keine Scherspannungen. Dies bedeutet, daß die Ablei-

tung der tangentialen Geschwindigkeit uτsenkrecht zur Wand (d.h. in Richtung der

Fl¨

achennormalen n) verschwindet.

un= 0 und ∂uτ

∂n = 0 .

3. Kontinuumsmechanische Modellierung

Die Schlupf-Randbedingung wird auch f¨

ur Symmetrieebenen angewendet. Im Falle ei-

ner Zweiphasenstr¨

omung kommt f¨

ur die Symmetrieebenen noch die Bedingung hinzu,

daß die Phasengrenze stets senkrecht auf eine Symmetrieebene auftrifft.

Einstr¨

ombedingung: Beide Geschwindigkeitskomponenten sind auf dem Rand fest

vorgegeben:

un=u0

nund uτ=u0

τ.

Ausstr¨

ombedingung: Beide Geschwindigkeitskomponenten ¨

andern sich nicht in der

Richtung senkrecht zum Rand:

∂un

∂n = 0 und ∂uτ

∂n = 0 .

4. Numerische Verfahren

Will man bei der Simulation zweiphasiger Str¨

omungen den Einfluß der Oberfl¨

achen-

spannung ber¨

ucksichtigen, so muß zu jedem Zeitpunkt die Form der Phasengrenze

bekannt sein, denn die Kr¨

ummung der Phasengrenze geht in die Impulsbilanz ein. Es

gibt verschiedene numerische Methoden, die Bewegung der Phasengrenze im Rechen-

gebiet zu beschreiben. Die Phasengrenze kann durch Punkte markiert werden, die

mit dem Geschwindigkeitsfeld transportiert werden. Ist die Phasengrenze dicht genug

mit Markierungspunkten belegt, so l¨

aßt sich aus der Position dieser Punkte stets die

Grenzfl¨

ache rekonstruieren. Eine andere Methode besteht darin, ein mitbewegtes Git-

ter zu verwenden, welches die Phasengrenze repr¨

asentiert. Im hier verwendeten Code

FS3D wird die Volume-of-Fluid (VoF) Methode verwendet.

4.1. Die VoF-Methode

Das VoF-Modell wurde von Hirt und Nichols (1981) vorgestellt. Seither wurde es

kontinuierlich verbessert (Rider und Kothe, 1998). Der VoF-Code FS3D wurde von

Rieber (2004) am ITLR Stuttgart entwickelt und verwendet die weiter unten beschrie-

bene PLIC-Methode zur Rekonstruktion der Phasengrenze.

Es wird von zwei inkompressiblen Fluiden ausgegangen, die sich auf molekularer Ebe-

ne nicht vermischen (l¨

osen). Dann ist der von beiden Fluiden eingenommene Gesamt-

raum Ω die disjunkte Vereinigung der Phasengrenze Γ mit den beiden Teilr¨

aume Ω+

und Ω−, wobei Ω+den von der dispersen Phase und Ω−den von der kontunierlichen

Phase eingenommenen Teilraum bezeichnet. Auf dem Gesamtraum Ω kann nun eine

sogenannte Phasenindikatorfunktion χdefiniert werden als

χ(x) = 









1 f¨

ur x∈Ω+

0 sonst .

(4.1)

Um das Str¨

omungsproblem numerisch l¨

osen zu k¨

onnen, muß der Raum mittels ei-

nes Gitters diskretisiert werden. Dies geschieht auch im hier angewendeten Finite-

4. Numerische Verfahren

Volumen-Verfahren, in welchem die integralen Bilanzgleichungen f¨

ur die Volumina

der einzelen Gitterzellen numerisch gel¨

ost werden. Aus der Phasenindikatorfunktion

kann auf dem Gitter in nat¨

urlicher Weise die sogenannte Volume-of-Fluid Variable f

definiert werden, die f¨

ur jede Zelle einen diskreten Wert annimmt. Die Verteilung der

Phasen wird implizit durch die sogenannte Volume-of-Fluid Variable fbeschrieben.

Jeder Zelle wird ein Wert fzugeordnet, der als Phasenindikator fungiert:

•f= 0 die Zelle enth¨

alt nur die kontinuierliche Phase .

•f= 1 die Zelle enth¨

alt nur die disperse Phase.

•0< f < 1 die Zelle enth¨

alt beide Phasen, mit dem Volumenanteil fan disperser

Phase.

Eine Zelle mit beiden Phasen (0 < f < 1) enth¨

alt demnach einen Teil der Phasen-

grenze. F¨

ur jeden Zeitschritt wird neben den Bilanzgleichungen f¨

ur Masse und Impuls

der Transport der Volumenanteilsfunktion fberechnet. Die Transportgleichung f¨

ur f

lautet:

∂f

∂t +∇·(fu) = 0 .(4.2)

Es gibt verschiedene M¨

oglichkeiten, aus einer gegebenen f-Verteilung die Phasen-

grenze zu rekonstruieren. Bei der “Piecewise Linear Interface Construction“ (PLIC)

Methode nach Puckett und Saltzmann (1992) wird innerhalb einer jeden Zelle ein ebe-

nes Fl¨

achenst¨

uck als N¨

aherung f¨

ur die Phasengrenze bestimmt. Dieses Fl¨

achenst¨

uck

ist so orientiert, daß die Fl¨

achennormale in Richtung ∇fzeigt. Zudem liegt das

Fl¨

achenst¨

uck so innerhalb der Zelle, daß sie die Zelle entsprechend dem f-Wert in

zwei Volumina f¨

ur die disperse und die kontinuierliche Phase teilt. Die Fl¨

achenst¨

ucke

benachbarter Zellen gehen nach dieser Vorgehensweise im Allgemeinen nicht stetig

ineinander ¨

uber. Die Kenntnis der Lage und Orientierung der Fl¨

achenst¨

ucke ist f¨

die Berechnung des konvektiven Transports der beiden Phasen erforderlich. Dies soll

an der Zelle mit dem VoF-Wert f= 0,15 gezeigt werden. Wenn in einem Zeitschritt

der Stoff um eine halbe Zellbreite nach oben transportiert wird, so gelangt aus dieser

Zelle nur kontinuierliche Phase (weiß) in die obere Nachbarzelle. Findet der Stoff-

transport stattdessen um eine halbe Zellbreite nach unten statt, so werden aus der

Zelle mit f= 0,15 Anteile beider Phasen in die untere Nachbarzelle transportiert.

4.2. Implementierung der Grenzfl¨

achenkraft

0,6

1 0,6

0,751

0,15 0

Abbildung 4.1.: Verteilung der VoF-Variablen in einem Gitter und Rekonstruktion der

Phasengrenze nach dem PLIC-Verfahren

4.2. Implementierung der Grenzfl¨

achenkraft

In die Grenzfl¨

achenkraft fΓgeht die Kr¨

ummung κΓein. Im Falle konstanter Ober-

fl¨

achenspannung gilt

fΓ=σ κΓnΓmit κΓ=−∇Γ·nΓ.(4.3)

Die Normale der Phasengrenze kann aus der Verteilung der VoF-Variablen bestimmt

werden:

nΓ=1

|∇f|∇f . (4.4)

Die Grenzfl¨

achenkraft fΓin (4.3) ist nur auf Γ definiert. Im hier verwendeten Code

werden alle Bilanzgleichungen in integraler Form gel¨

ost, wobei die betrachteten Vo-

lumina V, jeweils die Volumina der Zellen sind. Es bietet sich daher an, auch die

Grenzfl¨

achenkraft fΓin eine Volumenkraft fΓ,V zu ¨

uberf¨

uhren, sodaß

ZV∩Γ(t)

fΓdA =ZV

fΓ,V dV (4.5)

gilt. Man kann zeigen, daß im Falle konstanter Oberfl¨

achenspannung

fΓ,V =σ κΓ∇f(4.6)

die Bedingung (4.5) erf¨

ullt. Dies gilt, weil im Grenzfall unendlich kleiner Zellen

(V→0) der Ausdruck |∇f|in eine Dirac’sche Delta-Funktion ¨

ubergeht, die nur

an der Phasengrenze von Null verschieden ist. Die Berechnung der Grenzfl¨

achenkraft

4. Numerische Verfahren

in dieser Form hat jedoch den Nachteil, daß sie nicht konservativ ist. Das heißt, daß

die Impulserhaltung im numerischen Verfahren nicht gew¨

ahrleistet ist.

F¨

ur die in dieser Arbeit durchgef¨

uhrten Simulationen wurde die Grenzfl¨

achenkraft in

konservativer Form verwendet, d.h. sie ist als Divergenz eines Audrucks formuliert.

Im Falle konstanter Grenzfl¨

achenspannung gilt

∇Γ·(σI−σnΓ⊗nΓ) = σ κΓnΓ.(4.7)

Durch die Verwendung der auf (4.7) basierenden konservativen Form wird gew¨

ahrlei-

stet, daß auch in der numerischen Umsetzung die Impulserhaltung gilt. Eine genauere

Beschreibung dieser Methode findet sich bei (Lafaurie et al., 1994).

5. Dimensionsanalyse der Navier-Stokes-Gleichungen

Um die charakteristischen dimensionslosen Kennzahlen f¨

ur ein str¨

omungsmechani-

sches Problem zu ermitteln, k¨

onnen zwei unterschiedliche Methoden gew¨

ahlt werden.

Entweder man wendet das Π-Theorem an oder man formt die beschreibenden Glei-

chungen (hier: die Navier-Stokes-Gleichungen) so um, daß sie nur noch dimensionslose

Variablen enthalten. Im Folgenden soll der zweite Weg beschritten werden. In einem

ersten Schritt soll dieses Verfahren an den einphasigen Navier-Stokes-Gleichungen

demonstriert werden, um es danach auf die zweiphasigen Navier-Stokes-Gleichungen

auszudehnen.

5.1. Die einphasigen Navier-Stokes-Gleichungen in

dimensionsloser Form

Zun¨

achst sollen also die Navier-Stokes-Gleichungen f¨

ur einphasige Str¨

omungen unter-

sucht werden. Diese lauten f¨

ur inkompressible Medien

∂tu+u·∇xu=−1

ρ∇xp+µ

ρ∆xu+g,und (5.1)

∇x·u= 0 .(5.2)

Um besonders hervorzuheben, daß sich die r¨

aumlichen Ableitungen auf dimensions-

behaftete Ortskoordinaten beziehen, wird der Index ’x’ eingef¨

uhrt: Seien die Kom-

ponenten eines Ortsvektors x= (x1, x2, x3), so ist ∇x:= (∂/∂x1, ∂/∂x2, ∂/∂x3).Alle

Gr¨

oßen in (5.1) und (5.2) sind dimensionsbehaftet. Um die dimensionsbehafteten

Navier-Stokes-Gleichungen in dimensionslose Form zu bringen, werden die dimensi-

onsbehafteten Gr¨

oßen durch Normierung mit charakteristischen Gr¨

oßen in dimensi-

onslose Form ¨

uberf¨

uhrt. Beispielsweise kann die Geschwindigkeit uunter Verwendung

einer charakteristischen Geschwindigkeit Uin eine dimensionslose Geschwindigkeit

umgewandelt werden: v:= u/U. Charakteristische Gr¨

oßen sind in der Praxis so zu

w¨

ahlen, daß ihre Gr¨

oßenordung typisch f¨

ur die betreffende Anordnung ist. In einem

5. Dimensionsanalyse der Navier-Stokes-Gleichungen

Rohr beispielsweise bietet es sich an, als charakteristische Geschwindigkeit die mitt-

lere Geschwindigkeit U=˙

V /A zu w¨

ahlen, wobei Ader Rohrquerschnitt ist. Ebenso

k¨

onnen durch Verwendung einer charakteristischen Zeit τund einer charakteristischen

L¨

ange Ljeweils eine dimensionslose Zeit- und eine dimensionslose Ortskoordinate ein-

gef¨

uhrt werden. Es ergeben sich folgende drei dimensionslose Gr¨

oßen:

Zeit s:= t/τ ,

Ortskoordinate y:= x/L ,

Geschwindigkeit v:= u/U . (5.3)

Demnach bezeichnen hier vund ynicht die zweite Komponente der dimensionsbe-

hafteten Geschwindigkeit oder Ortskoordinate. Nun sollen die Terme in den dimensi-

onsbehafteten Navier-Stokes-Gleichungen sukzessive durch solche ersetzt werden, in

denen nur die dimensionslosen Gr¨

oßen vorkommen. Betrachten wir die Zeitableitung

der Geschwindigkeit:

v(s, y) = 1

Uu(τs, Ly)⇒∂v

∂s =τ

∂u

∂t ⇔∂u

∂t =U

∂v

∂s .(5.4)

Analog hierzu kann man die anderen Ableitungen behandeln. Beispielsweise ergibt

sich f¨

ur die Ortsableitungen des Druckes

∂p

∂xi

∂p

∂yi⇒ ∇xp=1

L∇yp . (5.5)

Insgesamt ergeben sich folgende Ersetzungen:

∂tu=U

τ∂sv;∇xu=U

L∇yv; (5.6)

∇xp=1

L∇yp′; ∆xu=U

L2∆yv.(5.7)

Dabei gelte p′(y) := p(x/L) . Die Gleichungen (5.6) und (5.7) werden nun in die

dimensionsbehafteten Navier-Stokes-Gleichungen (5.1) und (5.2) eingesetzt. Es folgt:

τ∂sv+U2

Lv·∇yv=−1

ρL ∇yp+µ

L2∆yv+g; (5.8)

L∇y·u= 0 .(5.9)

5.1. Die einphasigen Navier-Stokes-Gleichungen in dimensionsloser Form

Durch Multiplikation mit L/U kann Gl. (5.9) weiter vereinfacht werden zu ∇y·u= 0.

Es liegt nahe, auch die Vorfaktoren auf der linken Seite von Gl. (5.8) auf ¨

ahnliche

Weise zu beseitigen. Dies wird m¨

oglich, wenn die charakteristische Zeit τals als Quo-

tient aus charakteristischer L¨

ange und Geschwindigkeit definiert wird: τ:= L/U.

Eine solche Festlegung ist durchaus sinnvoll: Wenn beispielsweise Ldie L¨

ange einer

durchstr¨

omten Apparatur ist, und Ueine mittlere Str¨

omungsgeschwindigkeit, so er-

gibt sich aus L/U eine mittlere Aufenthaltszeit f¨

ur fluide Elemente in der Apparatur.

Mit dieser Festlegung kann auch die dimensionslose Impulsgleichung (5.8) vereinfacht

werden. Somit erh¨

alt man f¨

ur die dimensionslose Impulsgleichung

∂sv+v·∇yv=−1

ρ U2∇yp′+µ

ρ L U ∆yv+L

U2g.(5.10)

Bis auf den Druck p′sind alle Unbekannten in dieser Gleichung dimensionslos. Mit

der Einf¨

uhrung eines charakteristischen Druckes p0kann ein dimensionsloser Druck

˜p:= p′/p0eingef¨

uhrt werden. Die Vorfaktoren der Terme auf der rechten Seite werden

zu dimensionslosen Gr¨

oßen zusammengefaßt, die jeweils als Quotient aus zwei Kr¨

aften

interpretiert werden k¨

onnen.

Eulerzahl: Eu = p0/ρ U2Druck/Tr¨

agheit,

Reynoldszahl: Re = ρ U L/µ Tr¨

agheit/Viskosit¨

at,

Froudezahl: Fr = U2/gL Tr¨

agheit/Schwere.

Durch Einsetzen dieser Kennzahlen in die Impulsgleichung ergibt sich mit eg:= g/|g|

folgende Form f¨

ur die dimensionslose Navier-Stokes-Gleichungen:

∂sv+v·∇yv=−Eu ∇y˜p+1

Re ∆yv+1

Fr eg,(5.11)

∇y·v= 0 .(5.12)

Bei gegebenen Anfangs- und Randbedingungen ist die L¨

osung der einphasigen Navier-

Stokes-Gleichungen durch die drei Kennzahlen Re, Eu und Fr festgelegt. Mit anderen

Worten: Sind zwei durchstr¨

omte Apparaturen geometrisch ¨

ahnlich (d.h. bis auf Maß-

stabsvergr¨

oßerung formgleich), und sind Re, Eu und Fr identisch, so entwickeln sich

in ihnen auch ¨

ahnliche Str¨

omungsfelder. Zwei Str¨

omungsfelder sind ¨

ahnlich, wenn

es lineare Transformationen nach Art der Gleichungen (5.3) gibt, mit denen das ei-

ne in das andere ¨

uberf¨

uhrt werden kann. Dieser Zusammenhang wird benutzt, um

Versuche an maßstabsgetreuen Modellen im Str¨

omungskanal auf die Wirklichkeit zu

ubertragen.

5. Dimensionsanalyse der Navier-Stokes-Gleichungen

5.2. Die mehrphasigen Navier-Stokes-Gleichungen in

dimensionsloser Form

Analog zur Vorgehensweise im vorangegangenen Kapitel beginnen wir mit den di-

mensionsbehafteten Navier-Stokes-Gleichungen.

∂tu+u·∇xu=−1

ρ∇xp+µ

ρ∆xu+g,(5.13)

∇x·u= 0 ,(5.14)

[u] = 0 ,(5.15)

[pI−S]·n=σκ n.(5.16)

Gegen¨

uber den einphasigen Navier-Stokes-Gleichungen sind die Sprungbedingungen

(5.15) und (5.16) hinzugekommen. Durch (5.15) wird die Stetigkeit des Geschwindig-

keitsfeldes ¨

uber die Phasengrenze gefordert. Der Drucksprung an der Phasengrenze

durch die Oberfl¨

achenspannung wird durch (5.16) beschrieben. Die Gr¨

oßen Iund S

sind 3 ×3 Tensoren, wobei Ieine Einheitsmatrix ist. Sist der Spannungstensor:

S=µ(∇xu+∇xuT).(5.17)

Die Fl¨

achennormale der Phasengrenze nist eine Funktion des Ortes und der Zeit,

was durch die Schreibweise n(t, x) ausgedr¨

uckt werden soll. Wird zu dimensionslosen

Koordinaten s:= t/τ und y:= x/L ¨

ubergegangen, so ergibt sich eine neue Funktion

f¨

ur die Fl¨

achennormale:

n(s, y) := n(τs, Ly).(5.18)

Analog zum Druck in Gleichung (5.7) ergibt sich f¨

ur die Kr¨

ummung:

κ=−∇x·n=−1

L∇y·˜

n.(5.19)

Die Normalenvektoren nund ˜

nsind per definitionem dimensionslos. Die Kr¨

ummung

κentsteht durch eine Ortsableitung von nund hat dadurch die Dimension [1/L¨

ange].

Eine dimensionslose Kr¨

ummung kann also durch Multiplikation mit der charakteri-

stischen L¨

ange gewonnen werden:

˜κ:= L κ =−∇y·˜

5.2. Die mehrphasigen Navier-Stokes-Gleichungen in dimensionsloser Form

Die rechte Seite von Gleichung (5.16) kann also durch den Term σ˜κ˜

n/L ersetzt wer-

den, der von den dimensionslosen Koordinaten yund τabh¨

angt. Nun soll auch die

linke Seite als Funktion von dimensionslosen Koordinaten ausgedr¨

uckt werden. Unter

Verwendung von Gleichung (5.6) und p′:= p(τs, Ly) erh¨

alt die Sprungbedingung

folgende Form:

p′I−2µU

L˜

D·˜

n=1

Lσ˜κ˜

nmit ˜

D=1

2(∇yv+∇yvT).

Dies ist eine Gleichung in dimensionslosen Variablen, jedoch sind die Terme beider-

seits des Gleichheitszeichens noch dimensionsbehaftet. Durch Einf¨

uhrung eines cha-

rakteristischen Druckes p0kann ein dimensionsloser Druck ˜p:= p′/p0definiert werden.

Die Wahl eines beliebigen positiven ρ0und die Division der Sprungbedingung durch

ρ0U2ergibt eine dimensionslose Gleichung:

p0

ρ0U2˜pI−µ

ρ0UL 2˜

D˜

n=∆ρ

ρ0

∆ρ U2L˜κ˜

nmit ∆ρ=|ρc−ρd|.

Wird ein positives µ0gew¨

ahlt, so definiert ˜µ=µ/µ0eine dimensionslose Viskosit¨

und ˜

S= ˜µ˜

Deinen dimensionslosen Spannungstensor. Hiermit l¨

aßt sich die dimensi-

onslose Sprungbedingung umschreiben zu

Eu0˜pI−1

Re0

S˜

n=∆ρ

E¨

o˜κ˜

mit

Eu0=p0

ρ0U2; Re0=ρ0UL

µ0

; E¨

o = ∆ρU2L

σ.

Die dimensionslose Impulsgleichung (5.10) gilt auch im zweiphasigen Fall, jedoch sind

hier die Stoffeigenschaften ρund µdavon abh¨

angig, welche Phase am jeweiligen Ort

vorliegt. Durch Multiplikation von (5.10) mit ρ/ρ0ergibt sich:

ρ0

∂sv+ρ

ρ0

v·∇yv=−1

ρ U2∇yp′+µ

ρ0L U ∆yv+ρ

ρ0

U2g.

Mit ˜ρ=ρ/ρ0, ˜p=p/p0und ˜µ=µ/µ0folgt:

˜ρ∂sv+ ˜ρv·∇yv=−p0

ρ0U2∇y˜p+ ˜µµ0

ρ0L U ∆yv+ ˜ρL

U2g;

⇒˜ρ∂sv+ ˜ρv·∇yv=−Eu0∇y˜p+ ˜µ1

Re0

∆yv+ ˜ρ1

Freg.

5. Dimensionsanalyse der Navier-Stokes-Gleichungen

Handelt es sich um eine Zweiphasenstr¨

omung mit einer kontinuierlichen und einer

dispersen Phase, so werden in der Regel die Stoffeigenschaften der kontinuierlichen

Phase ρcund µcals Bezugsgr¨

oßen verwendet: ρ0:= ρc;µ0:= µc. Im Folgenden werden

mit Eu und Re die auf die kontinuierliche Phase bezogenen Kennzahlen bezeichnet.

Damit ergeben sich die dimensionslosen zweiphasigen Navier-Stokes-Gleichungen:

˜ρ∂sv+ ˜ρv·∇yv=−Eu ∇y˜p+ ˜µ1

Re ∆yv+ ˜ρ1

Fr eg; (5.20)

∇yv= 0 ; (5.21)

[v] = 0 ; (5.22)

Eu ˜pI−1

Re ˜

S˜

n=∆ρ

E¨

o˜κ˜

n; (5.23)

mit

Eu = p0

ρcU2; Re = ρcUL

µ0

; E¨

o = ∆ρU2L

σ; Fr = U2

g L .

Diese vier Konstanten sind dimensionslos. Die Sprungbedingung enth¨

alt zus¨

atzlich die

beiden dimensionslosen Gr¨

oßen ∆ρ/ρcund das Viskosit¨

atenverh¨

altnis µd/µc, welches

als Faktor im Spannungstensor ˜

Senthalten ist. Durch die Wahl von p0:= ρcU2wird

Eu zu Eins. Durch die Wahl von U= (g L)1/2wird Fr zu Eins. F¨

ur E¨

o ergibt sich

damit:

E¨

o = g∆ρ L2

σ.

In dieser Gestalt heißt E¨

o auch ”E¨

otv¨

os-Zahl“ und kann als Quotient von Auftriebs-

kraft und Oberfl¨

achenspannungskraft eines fluiden Partikels mit der charakteristi-

schen L¨

ange Linterpretiert werden. Es verbleiben vier dimensionslose Kennzahlen

zur Charakterisierung zweiphasiger Str¨

omungen:

∆ρ

ρc

;λ=µd

µc

; Re = ρcUL

µ0

; E¨

o = g∆ρ L2

σ.(5.24)

Anstelle von ∆ρ/ρcwird oftmals die Morton-Zahl verwendet, welche die Stoffeigen-

schaften der beiden Phasen charakterisiert:

Mo = g∆ρ µ4

ρ2

cσ3.

Wegen

∆ρ

ρc2

=E¨

Mo Re4

sind ∆ρ/ρcund Mo im Kennzahlensatz (5.24) austauschbar.

6. Validierung des Codes FS3D

6.1. Die Hadamard-Rybczinski-L¨

osung

Eine analytische L¨

osung der Navier-Stokes-Gleichungen f¨

ur ein Zweiphasenproblem

ist im Allgemeinen nicht m¨

oglich. Im Grenzfall Re →0 (d.h. bei kleinen Geschwindig-

keiten u,”Kriechstr¨

omung“) kann der nichtlineare Term u∇uin der Impulsgleichung

(5.13) vernachl¨

assigt werden, weil hier uquadratisch auftritt. Es ergibt sich eine ver-

einfachte Impulsgleichung f¨

ur Kriechstr¨

omungen, die Stokes-Gleichung genannt wird:

∂tu=−1

ρ∇xp+µ

ρ∆xu+g.(6.1)

Mit dieser Vereinfachung existiert eine analytische L¨

osung f¨

ur fluide Kugeln, die von

Hadamard und Rybzinski gefunden wurde. F¨

ur das Geschwindigkeitsfeld lautet die

L¨

osung in Kugelkoordinaten außerhalb der fluiden Kugel:

ur(r, θ) = −wT"1−R

3λ+ 2λ

2λ+ 2 +R

r2λ

2λ+ 2#cos θ; (6.2)

uθ(r, θ) = wT"1−R

3λ+ 2λ

4λ+ 4 −R

r3λ

λ+ 1#sin θ; (6.3)

und innerhalb der fluiden Kugel:

ur=wT

2λ+ 2 1−r

R2cos θ; (6.4)

uθ=−wT

2λ+ 2 1−2r

R2sin θ(6.5)

mit dem Kugelradius Rund dem Viskosit¨

atenverh¨

altnis λ. Dabei sind urund uθ

jeweils die radiale und die vertikale Relativgeschwindigkeit bezogen auf die Schwer-

punktbewegung, die mit der Geschwindigkeit wTin Richtung θ= 0 erfolgt. Die

L¨

osung ist axialsymmetrisch mit der Linie θ= 0 bzw. θ=πals Symmetrieachse. Die

Ebene θ=π/2 - im Folgenden ”¨

Aquatorebene“ genannt - verl¨

auft durch den Mittel-

punkt und senkrecht zur Bewegungsrichtung des fluiden Partikels. Die ¨

Aquatorebene

6. Validierung des Codes FS3D

ist dadurch ausgezeichnet, daß auf ihr die radiale Geschwindigkeitskomponente ur

sowohl innerhalb als auch außerhalb der fluiden Kugel verschwindet. Das Geschwin-

digkeitsfeld innerhalb der fluiden Kugel ist ein torusf¨

ormiger Wirbelring, der auch

Hill’scher Wirbel genannt wird. Im Wirbelzentrum gilt ur=uθ= 0. Diese Bedin-

gung ist auf der ¨

Aquatorebene bei r=R/√2 erf¨

ullt. Weitere ausgezeichnete Punkte

auf der ¨

Aquatorebene sind das Kugelzentrum (r= 0) und der ¨

Aquator (r=R). Im

Kugelzentrum ist die Geschwindigeit ucsenkrecht zur ¨

Aquatorebene:

uc=wT

2λ+ 2 .(6.6)

Auf dem ¨

Aquator ist die Geschwindigkeit senkrecht zur ¨

Aquatorebene −uc. Insge-

samt ist das Geschwindigkeitsprofil innerhalb der fluiden Kugel auf der ¨

Aquatorebene

gem¨

aß (6.5) parabolisch. Ist das Geschwindigkeitsfeld bekannt, kann das Druckfeld

durch Einsetzen der Geschwindigkeit in die Navier-Stokes-Gleichungen ermittelt wer-

den. Im Anhang A.2 wird dies am Beispiel des Hill’schen Wirbels gezeigt. Insgesamt

ergibt sich f¨

ur die fluide Kugel folgendes Druckfeld:

p=p0−ρcg r cos θ+wTµc

RR

r23λ+ 2

2λ+ 2 cos θf¨

ur r > R ; (6.7)

p=p0−ρdg r cos θ−5wTµc

λ+ 1 cos θ+2σ

Rf¨

ur r < R . (6.8)

Durch Integration des durch (6.7) gegebenen Außendruckes entlang des Kugelrandes

(r=R) erh¨

alt man den Druckbeitrag zum Widerstandsbeiwert

CD=8

3 Re 2 + 3λ

1 + λ.(6.9)

Durch Integration von S·erentlang des Kugelrandes ergibt sich der Beitrag viskoser

Kr¨

afte zum Widerstandsbeiwert

CD,S =16

3 Re 2 + 3λ

1 + λ.(6.10)

Die Summe dieser Beitr¨

age ergibt den gesamten Widerstandsbeiwert auf eine fluide

Kugel zu

CD,p =8

Re 2 + 3λ

1 + λ.(6.11)

6.2. Numerische L¨

osungen f¨

ur langsame, sph¨

arische fluide Partikel

F¨

ur eine Gasblase gilt λ≈0 somit CD= 16/Re. Eine feste Kugel kann als eine fluide

Kugel mit unendlich hoher Viskosit¨

at angesehen werden. Durch den Grenz¨

ubergang

λ→ ∞ ergibt sich der Widerstandsbeiwert f¨

ur feste Kugeln als CD= 24/Re. Beide

Beziehungen gelten gem¨

aß Voraussetzung nur f¨

ur Re ≪1. Durch Gleichsetzen der

Widerstandskraft mit der Auftriebskraft kann aus (12.4) die Aufstiegsgeschwindigkeit

ermittelt werden. Diese lautet:

wT=2

gR2(ρc−ρd)

µc

1 + λ

2 + 3λ.(6.12)

6.2. Numerische L¨

osungen f¨

ur langsame, sph¨

arische fluide

Partikel

Die L¨

osung von Hadamard-Rybczinski f¨

ur langsam aufsteigende/sinkende fluide Par-

tikel (d.h. Re≪1) kann zum Vergleich mit numerischen L¨

osungen herangezogen wer-

den, um die Zuverl¨

assigkeit des verwendeten Codes zu ¨

uberpr¨

ufen. Man spricht in die-

sem Zusammenhang von der ”Validierung“ eines Codes. Um den Code FS3D zu vali-

dieren, wird zun¨

achst ein Modellsystem gew¨

ahlt, welches die Simulation eines langsam

aufsteigenden fluiden Partikels erm¨

oglicht. Durch eine hohe Viskosit¨

at µcder kontinu-

ierlichen Phase kann die Reynoldszahl klein gehalten werden. F¨

ur das Modellsystem

wurden folgende Stoffdaten gew¨

ahlt: µc=0,2 Pa s, ρc= 2000 kg/m3,ρd= 1000 kg/m3

und σ= 72 mN/m. Die Oberfl¨

achenspannung σist identisch mit der von Wasser an

Luft. F¨

ur die disperse Phase werden verschiedene Viskosit¨

aten µdverwendet so daß

das Viskosit¨

atenverh¨

altnis λ=µd/µcim Intervall 0,125 ≤λ≤2 liegt. F¨

ur ver-

schiedene λwird die Endgeschwindigkeit eines Tropfens mit einem Durchmesser von

d=0,2 cm berechnet. Das Rechengebiet ist quaderf¨

ormig mit vertikalen W¨

anden im

Abstand von 8dbzw. 16d, d.h. es wird der Aufstieg eines Tropfens in einem Rohr mit

einem quadratischen Querschnitt simuliert. Die Kantenl¨

ange des Rohrquerschnitts

betr¨

agt entweder 8doder 16d. Die Aufstiegsgeschwindigkeit sinkt mit wachsendem λ,

d.h., mit der Erh¨

ohung von µc(siehe Abbildung 6.1). Grunds¨

atzlich sind die nu-

merisch ermittelten Aufstiegsgeschwindigkeiten kleiner als nach Gleichung (6.12).

Diese Abweichung wird mit großer Wahrscheinlichkeit durch die Seitenw¨

ande ver-

ursacht, denn mit gr¨

oßerem Wandabstand dWsteigen die Endgeschwindigkeiten wT.

Im Grenzfall dW→ ∞ ist das kontinuierliche Fluid unendlich ausgedehnt, was Vor-

aussetzung f¨

ur die L¨

osung von Hadamard-Rybzinski ist. Aufstiegsgeschwindigkeiten

6. Validierung des Codes FS3D

von 1,0< wT<1,4 entsprechen 0,1<Re <0,14, womit die zweite Voraussetzung

Re ≪1 erf¨

ullt wird. Um das numerische Ergebnis f¨

ur das Geschwindigkeitsfeld im Fall

λ= 1 mit der analytischen L¨

osung zu vergleichen, wird die senkrechte Komponente

der Relativgeschwindigkeit wr=w−wTauf einer horizontalen Linie in x-Richtung

durch den Blasenschwerpunkt untersucht. Diese Linie liegt auf der ¨

Aquatorebene,

wo horizontale Geschwindigkeitskomponenten Null sind. Zun¨

achst wird die analyti-

sche L¨

osung in kartesische Koordinaten umgewandelt: Sei der Blasenschwerpunkt bei

x=0, so gilt auf der ¨

Aquatorebene: wr(x) = uθ(r=x, θ =π/2). Hiermit k¨

onnen

die Werte f¨

ur wrleicht aus (6.5) und (6.3) hergeleitet werden. Wie aus Abbildung

6.2 hervorgeht, ist der Nulldurchgang wr= 0 f¨

ur die theoretische L¨

osung wie f¨

ur die

analytische L¨

osung an derselben Stelle. Hier liegt das Zentrum des Hill’schen Wirbels

(vgl. Seite 28). In der analytischen L¨

osung gilt wr(x→ ∞) = −wT, denn in unendlich

großer Entfernung vom Tropfen ruht die kontinuierliche Phase (w= 0 ⇒wr=−wT).

In den Simulationen hingegen ist das Rechengebiet durch vertikale W¨

ande begrenzt.

Die durch den aufsteigenden Tropfen verdr¨

angte kontinuierliche Phase muß zwischen

dem Tropfenrand und den W¨

anden abw¨

arts str¨

omen. Daher sind außerhalb des Trop-

0,0 0,5 1,0 1,5 2,0

0,0

0,2

0,4

0,6

0,8

1,0

1,2

1,4

1,6

[cm/ s]

=16d

=8d

Abbildung 6.1.: Endgeschwindigkeit simulierter Modelltropfen in Abh¨

angigkeit vom Vis-

kosit¨

atenverh¨

altnis λ. Kurve: wTnach der L¨

osung von Hadamard-Rybczinski, Gleichung

(6.12). Gestrichelte Linie: wTf¨

ur λ→ ∞ nach Gleichung (6.12). Symbole: Simulierte Auf-

stiegsgeschwindigkeiten wT; Kreise: Wandabstand dW= 8d. Dreiecke: dW= 16d.

6.2. Numerische L¨

osungen f¨

ur langsame, sph¨

arische fluide Partikel

0 1 2 3 4 5 6 7 8

-1,2

-1,0

-0,8

-0,6

-0,4

-0,2

0,0

0,2

0,4

=16d

= 8d

/ w

x / R

Abbildung 6.2.: Senkrechte Komponente der Relativgeschwindigkeit wr(normiert auf

die Aufstiegsgeschwindigkeit wT) in Abh¨

angigkeit vom horizontalen Abstand x zum Trop-

fenschwerpunkt (normiert auf den Tropfenradius d/2) f¨

ur λ= 1. Kurve: Analytische

L¨

osung nach Hadamard-Rybczinski. Kreise: Wandabstand dW= 8d. Dreiecke: Wandab-

stand dW= 16d. Gestrichelte vertikale Linie: Lage der Phasengrenze.

6. Validierung des Codes FS3D

fens die Betr¨

age von wrf¨

ur die numerischen L¨

osungen gr¨

oßer als in der analytischen

L¨

osung. Auch innerhalb des Tropfens sind die numerischen Betr¨

age von wrgr¨

oßer.

Exemplarisch f¨

ur die St¨

arke der inneren Zirkulation wird die Geschwindigkeit im

Tropfenzentrum n¨

aher untersucht. F¨

ur die vertikale Relativgeschwindigkeit im Trop-

fenzentrum wr,c =wr(x= 0) folgt aus (6.6)

wr,c

1 + λ.(6.13)

F¨

ur das eingangs beschriebene Modellsystem wird wr,c numerisch ermittelt und mit

(6.13) verglichen (Abbildung 6.3). Das numerische Resultat ist stets gr¨

oßer als der

analytische Wert, jedoch sinkt der Unterschied mit zunehmendem Wandabstand.

Nachdem f¨

ur die Geschwindigkeitsfelder eine befriedigende ¨

Ubereinstimmung zwi-

schen numerischen und analytischen L¨

osungen gefundenen wurde, sollen nun die

Druckfelder verglichen werden. Auf der ¨

Aquatorebene (θ= 0) ist der Druck gem¨

aß

der Gleichungen (6.7) und (6.8) in beiden Fluiden konstant. Lediglich an der Phasen-

grenze tritt der durch die Oberfl¨

achenspannung verursachte Drucksprung 2σ/R auf.

In der numerischen L¨

osung ist der Drucksprung ¨

uber mehrere Zellen verschmiert, weil

die Oberfl¨

achenspannung als Volumenkraft modelliert ist, die auch in der Umgebung

0,0 0,5 1,0 1,5 2,0

0,0

0,1

0,2

0,3

0,4

0,5

0,6

1/(2+2 )

= 8d

=16d

c,r

/ w

Abbildung 6.3.: Senkrechte Komponente der Relativgeschwindigkeit im Tropfenschwer-

punkt in Abh¨

angigkeit vom Viskosit¨

atenverh¨

altnis λ. Kurve: Analytische L¨

osung nach

Hadamard-Rybczinski. Kreise: Wandabstand dW= 8d. Dreiecke: Wandabstand dW= 16d.

6.2. Numerische L¨

osungen f¨

ur langsame, sph¨

arische fluide Partikel

der Phasengrenze pr¨

asent ist (Abbildung 6.4). Als weitere Kontroll-Linie bietet sich

die Vertikale durch den Tropfenschwerpunkt an. Gem¨

aß (6.8) ¨

andert sich hier der

Innendruck des Tropfens linear mit der vertikalen Koordinate z=rcos θ:

p=p0−ρdg z −5wTµc

λ+ 1 z+2σ

Rf¨

ur |z|< R . (6.14)

F¨

ur den Außendruck gilt auf der vertikalen Kontroll-Linie:

p=p0−ρcg z +wTµcR

z|z|

3λ+ 2

2λ+ 2 f¨

ur |z|> R . (6.15)

In Abbildung 6.5 ist der Gesamtdruck f¨

ur einen Tropfen entlang der Vertikalen ab-

gebildet. Der hydrodynamische Druckanteil ist jeweils der Gesamtdruck abz¨

uglich

der konstanten Beitr¨

age und des hydrostatischen Druckes −ρ g z. Der dynamische

Druckanteil pDaußerhalb des Tropfens entlang der Vertikalen ist in Abbildung 6.6

dargestellt. Der dynamische Druck pDinnerhalb des Tropfens hat die Steigung

∂

∂z pD=−5wTµc

λ+ 1 .(6.16)

0,0 0,5 1,0 1,5 2,0

100

120

140

160

p [Pa]

x / R

Abbildung 6.4.: Druck entlang einer horizontalen Linie durch den Tropfenschwerpunkt.

Die horizontale Verschiebung x ist durch den Tropfenradius R normiert (R=1 mm; λ= 1).

Punkte: Numerische L¨

osung. Linie: Analytische L¨

osung.

6. Validierung des Codes FS3D

-2 -1 0 1 2

-50

100

150

200

p [Pa]

z / R

Abbildung 6.5.: Druck in Abh¨

angigkeit von der H¨

ohe auf einer vertikalen Kontroll-Linie

durch den Tropfenmittelpunkt (d=2 mm, λ= 1).

-5 -4 -3 -2 -1 0 1 2 3 4 5

-3

-2

-1

p [Pa]

z / R

Abbildung 6.6.: Dynamischer Druck in Abh¨

angigkeit von der H¨

ohe auf einer vertikalen

Kontroll-Linie durch den Tropfenmittelpunkt (d=2 mm, λ= 1, Wandabstand dW= 8d).

6.2. Numerische L¨

osungen f¨

ur langsame, sph¨

arische fluide Partikel

Durch Einsetzen von (6.12) in (6.16) folgt:

∂

∂z pD=−10

3(ρc−ρd)λ

2 + 3λ.(6.17)

F¨

ur verschiedene Viskosit¨

atenverh¨

altnisse λwird die Steigung ∂pD/∂z im Tropfenin-

neren aus numerischen L¨

osungen ermittelt und mit dem analytischen Wert nach (6.17)

verglichen (Abbildung 6.7). Insgesamt wird das Druckfeld innerhalb und außerhalb

des Tropfens in der numerischen L¨

osung befriedigend wiedergegeben. Lediglich direkt

an der Phasengrenze kommt es aufgrund der Modellierung der Oberfl¨

achenspannung

als Volumenkraft zu gr¨

oßeren Abweichungen von der analytischen L¨

osung.

0,0 0,5 1,0 1,5 2,0

- d/dz (p

) [Pa / mm]

Abbildung 6.7.: Steigung des dynamischen Druckes im Tropfen (d=2 mm) in Abh¨

angig-

keit vom Viskosit¨

atenverh¨

alntis λ. Kurve: Analytische L¨

osung (Gleichung 6.16). Symbole:

Numerische L¨

osung. Kreise: Wandabstand dW= 8d. Dreiecke: Wandabstand dW= 16d.

6. Validierung des Codes FS3D

Teil I.

Bewegung von ¨

Oltropfen in Wasser

Motivation

Die Bewegung von ¨

Oltropfen in Wasser ist ein Problem, welches bei der F¨

orderung,

der Verarbeitung und dem Transport von Roh¨

ol eine große Rolle spielt. Durch die Ver-

brennung fossiler Energietr¨

ager wird Kohlendioxyd erzeugt, welches zur Erw¨

armung

der Erdatmosph¨

are wesentlich beitr¨

agt. In den letzten Jahren werden vermehrt An-

strengungen unternommen, Kohlendioxyd aus Abgasen abzuscheiden, um eine Frei-

setzung in der Atmosp¨

are zu verhindern. Dabei entsteht das Problem der Endlagerung

des Gases, welches in hohen Konzentrationen t¨

odlich ist. Eine zeitlang wurde erwogen,

Kohlendioxyd in der Tiefsee endzulagern. Unter den dort vorherrschenden Dr¨

ucken

ist es fl¨

ussig und dichter als Wasser, so daß es ab einer bestimmten Tiefe von selbst

bis zum Meeresboden absinken w¨

urde. Inzwischen bevorzugt man jedoch die L¨

osung,

Kohlendioxyd in ausgebeutete ¨

Olfelder zu pumpen. Ein großer Teil der ¨

Olfelder liegt

unter dem Meeresgrund. Daher ist es auch f¨

ur die Planung solcher Projekte wichtig,

das Verhalten fl¨

ussigen Kohlendioxyds in der Tiefsee zu kennen. Vor diesem Hinter-

grund wurden von Haljasmaa (2006) Experimente zur Sink- beziehungsweise Auf-

stiegsgeschwindigkeit von Tropfen fl¨

ussigen Kohlendioxyds und anderer Fl¨

ussigkeiten

in Wasser unternommen. In Wasser eignet sich ¨

Ol als Modellsubstanz f¨

ur fl¨

ussiges

Kohlendioxyd, weil die Stoffeigenschaften (wie z.B. die Morton-Zahl) beider Systeme

ahnlich sind. W¨

ahrend Kohlendioxyd nur unter hohem Druck fl¨

ussig bleibt, k¨

onnen

Oltropfen unter Normalbedingungen untersucht werden und sind daher leichter hand-

habbar.

Haljasmaa hat in seinen Experimenten Mais¨

ol und Rizinus¨

ol verwendet. Die Stoff-

eigenschaften der ¨

Olsorten und von Wasser sind in Tabelle 6.1 aufgelistet. F¨

ur die

Olsorten sind die Oberfl¨

achenspannungen σgegen¨

uber Wasser eingetragen. Zum Ver-

gleich: Die Oberfl¨

achenspannung f¨

ur Wasser gegen¨

uber Luft betr¨

agt 72 mN/m.

Fluid ρl[kg/m3]µl[mPa s] σ[mN/m] log10(Mo)

Mais¨

ol 917 45 20,1 9E-11

Rizinus¨

ol 957 900 16,8 7E-11

Wasser 998 1 - -

Tabelle 6.1.: Stoffdaten f¨

ur verschiedene Fluide aus Haljasmaa (2006)

7. Simulation von Mais¨

oltropfen

7.1. Konfiguration des Rechengebietes

Ziel der Rechnungen ist es, die Bewegung in einem ”unendlich“ ausgedehnten Was-

servolumen zu simulieren. Andererseits kann nur ein begrenztes Gebiet numerisch re-

pr¨

asentiert werden. Der hier verwendete Code FS3D sieht ausschließlich quaderf¨

ormi-

ge Rechengebiete vor. Die Randbedingungen sind so zu w¨

ahlen, daß die r¨

aumliche Be-

grenzung m¨

oglichst geringen Einfluß auf die Tropfenbewegung hat. So wird an den Sei-

tenw¨

anden des Rechengebietes eine Schlupf-Randbedingung verwendet, um Reibung

an der Wand auszuschließen. Am oberen Rand wird eine konstante Geschwindigkeit

w0vorgegeben, um den Tropfen durch Gegenstrom m¨

oglichst lange im Rechengebiet

zu halten. Am unteren Rand wurde eine Ausstr¨

om-Randbedingung gew¨

ahlt, wodurch

die oben eingef¨

uhrte Fl¨

ussigkeit das Rechengebiet wieder verlassen kann. Wenn nur

zweidimensionale oder eindimensionale Trajektorien zugelassen werden, so kann das

Rechengebiet durch Verwendung einer oder zweier Symmetrieebenen verkleinert wer-

den. Dadurch verringert sich der Rechenaufwand erheblich. Die Symmetrieebenen

verlaufen durch den Blasenschwerpunkt. F¨

ur zweidimensionale Trajektorien wird die

x-z-Ebene als Symmetrieebene verwendet, f¨

ur eindimensionale Trajektorien zus¨

atz-

lich die y-z-Ebene. Das Gitter ist kartesisch mit kubischen Zellen. Als Wandabstand

dWwird der Abstand zweier gegen¨

uberliegender W¨

ande mit Schlupf-Randbedingung

bezeichnet.

7.2. Einfluß des Wandabstandes auf die Aufstiegsgeschwindigkeit

Ein Mais¨

oltropfen mit einem ¨

Aquivalenzdurchmesser von d=5,7 mm steigt im Ex-

periment geradlinig auf. Daher k¨

onnen in der Simulation dieses Vorgangs ohne Ein-

schr¨

ankung der Physikalit¨

at zwei Symmetrieebenen verwendet werden. In Abbildung

7.1 ist die Endgeschwindigkeit wTdes ¨

Oltropfens in Abh¨

angigkeit vom Wandabstand

DWdargestellt. F¨

ur den Wandeinfluß auf die Aufstiegsgeschwindigkeit von Blasen

7. Simulation von Mais¨

oltropfen

und Tropfen im zylindrischen Rohr mit Durchmesser Dschlagen Clift et al. (1978)

folgende Korrelation vor:

wT,∞

="1−d

D2#3/2

.(7.1)

Hierbei ist wT,∞die Aufstiegsgeschwindigkeit ohne Wandeinfluß (D=∞). Es wird

E¨

o<40, Re >0 und d/D < 0,6 vorausgesetzt. Diese Bedingungen sind f¨

ur den

Mais¨

oltropfen mit d=5,7 mm erf¨

ullt. Setzt man D=DWein, so kann f¨

ur jeden Wand-

abstand aus dem jeweils berechneten wTein theoretischer Wert f¨

ur wT,∞ermittelt

werden. Dieser Wert stellt jedoch aus zwei Gr¨

unden lediglich eine obere Absch¨

atzung

dar: Erstens ist der Kanalquerschnitt in der Simulation quadratisch, und hat damit

eine gr¨

oßere Querschnittsfl¨

ache als ein Kreisrohr mit dem Durchmesser DW. Zwei-

tens wird durch die Schlupf-Randbedingung in der Simulation Reibung an der Wand

unterbunden. Wenn die Schlupfbedingung an den W¨

anden durch eine Haftbedingung

ersetzt wird, ¨

andert sich sich die Aufstiegsgeschwindigkeit der Tropfen in den Simu-

lationen allerdings nicht signifikant. Wegen des quadratischen Rohrquerschnitts wird

wT,∞durch Anwendung von (7.1) auf die simulierten Aufstiegsgeschwindigkeiten wT

systematisch ¨

ubersch¨

atzt. Der so ermittelte Wert kann jedoch als obere Schranke

interpretiert werden. In Abbildung 7.1 konvergieren wT,∞und wTschnell. Erh¨

oht

man den Wandabstand von 8dauf 16d, so ¨

andert sich die Aufstiegsgeschwindigkeit

kaum. Daher wird f¨

ur die folgenden Simulationen dW= 8dals Standardwandabstand

gew¨

ahlt.

7.3. Geschwindigkeiten innerhalb des Tropfens

Die Viskosit¨

at von ¨

Ol ist in der Regel um ein Vielfaches h¨

oher als die von Wasser. So

ist f¨

ur das hier untersuchte Mais¨

ol das Viskosit¨

atsverh¨

altnis λ= 45, f¨

ur das Castor¨

betr¨

agt es sogar 900. Im Grenzfall λ→ ∞kommen innere Geschwindigkeiten in einem

Tropfen zum Erliegen und er verh¨

alt sich wie ein festes Partikel. Bei großem, jedoch

endlichem λist zu erwarten, daß die Relativgeschwindigkeiten im fluiden Partikel

sehr klein sind. Dadurch herrscht an der Phasengrenze in tangentialer Richtung quasi

Haftbedingung. Jedoch ist eine langsame Deformation des Partikels durch geringe Ge-

schwindigkeiten der Phasengrenze in normaler Richtung auch bei großen λm¨

oglich.

7.3. Geschwindigkeiten innerhalb des Tropfens

T,inf

/ d

Abbildung 7.1.: Endgeschwindigkeit simulierter Mais¨

oltropfen (d=5,7 mm) in Abh¨

angig-

keit vom Wandabstand (Quadrate). Kreise: wT,∞nach (7.1)

Um eine Absch¨

atzung zu erhalten, ab welchem Viskosit¨

atenverh¨

altnis fluide Parti-

kel als ”fest“ im Sinne der tangentialen Geschwindigkeit an der Phasengrenze sind,

vergleicht Haljasmaa (2006) die Standardkurve f¨

ur den Widerstandsbeiwert fester

Kugeln mit einer Korrelation von Feng und Michaelides f¨

ur den Widerstandsbeiwert

sph¨

arischer fluider Partikel ohne Surfactants.

CD(Re, λ) = 









2−λ

2CD(Re,0) + 4λ

6 + λCD(Re,2) f¨

ur 0 ≤λ≤2 ;

λ+ 2 CD(Re,2) + λ−2

λ+ 2 CD(Re,∞) f¨

ur 2 ≤λ≤ ∞ ;

(7.2)

mit

CD(Re,0) = 48

Re 1 + 2,21

√Re −2.14

Re ,

CD(Re,2) = 17,0 Re−2/3,

CD(Re,∞) = 24

Re 1 + 1

6Re2/3.

Diese Korrelation f¨

ur den Bereich Re <1000 erhalten Feng und Michaelides mittels

einer numerischen Untersuchung, die eine Erweiterung der L¨

osung von Hadamard

und Rybczinski hin zu mittleren Reynoldszahlen darstellt, indem der Ansatz von

Levich (1962) auf λ > 0 verallgemeinert wird. In Abbildung 7.2 ist die Korrelation

7. Simulation von Mais¨

oltropfen

(7.2) f¨

ur verschiedene Werte von λsowie die Standard-Drag-Kurve f¨

ur feste Kugeln

nach Gleichung (2.6) eingezeichnet. Ab λ≥20 stellt die Korrelation im Bereich

50 <Re <1000 eine gute Approximation der Standard-Drag-Kurve dar. Weil eine

weitere Erh¨

ohung des Viskosit¨

atenverh¨

altnisses ¨

uber 20 hinaus zu keiner Erh¨

ohung

des Widerstandsbeiwertes f¨

uhrt, geht Haljasmaa davon aus, daß ab λ= 20 die Pha-

sengrenze in tangentialer Richtung unbeweglich ist. F¨

ur λ= 0 hingegen ergibt sich

aus (7.2) der Drag f¨

ur kugelf¨

ormige Blasen in reinem Wasser. Ab Re = 500 steigt der

Drag f¨

ur die Blasen stark an, was auf die zunehmende Deformation der Blasen und

das Einsetzen horizontaler Oszillationen zur¨

uckzuf¨

uhren ist. Diese Ph¨

anomene sind

in der Korrelation nicht ber¨

ucksichtigt. Nach der Annahme von Haljasmaa m¨

ußten

bei in Wasser aufsteigenden Mais¨

oltropfen (λ= 45) die Geschwindigkeiten innerhalb

des Tropfens und damit auch die tangentiale Geschwindigkeitskomponente an der

Phasengrenze vernachl¨

assigbar klein sein.

Um diese Vermutung zu ¨

uberpr¨

ufen, wird f¨

ur eine Serie von Tropfen mit gleichem

Durchmesser (d=5,7 mm) und verschiedenen Viskosit¨

aten µdder freie Aufstieg in

Wasser simuliert. F¨

ur die Oberfl¨

achenspannung und Dichte werden die Stoffwerte

Abbildung 7.2.: Dragkoeffizient in Abh¨

angigkeit von der Reynoldszahl f¨

ur verschiedene

Partikel: Fette Kurve: Standard-Drag f¨

ur feste Kugeln. D¨

unne Linien: Korrelation (7.2) f¨

verschiedene Viskosit¨

atenverh¨

altnisse λ. Offene Symbole: Experimentell bestimmter Drag

f¨

ur Luftblasen in reinem Wasser. A: Maxworthy et al. (1996) B: Haberman und Morton

(1953). Fette gestrichelte Linie: Schirmblasen, CD= 8/3.

7.3. Geschwindigkeiten innerhalb des Tropfens

0.00

0.10

0.20

0.30

0.40

0.50

0.60

Abbildung 7.3.: Geschwindigkeitsfeld eines Mais¨

oltropfens (λ= 45, d=5,7 mm) nach

Erreichen der Endgeschwindigkeit wT. (Bild um 90◦gedreht, Aufstiegsrichtung nach links.)

von Mais¨

ol ¨

ubernommen (σ= 20,1 mN/m, ρ= 917 kg/m3). Aus Abbildung 7.4

ist erkennbar, daß die Aufstiegsgeschwindigkeit auch ¨

uber λ= 20 mit steigendem λ

sinkt. F¨

ur Tropfen mit λ= 20 und λ= 45 wird das Geschwindigkeitsfeld im Tropfen-

inneren nach Erreichen der Endgeschwindigkeit wTn¨

aher untersucht. Betrachtet man

die vertikale Relativgeschwindigkeit wr=w−wTauf einer horizontalen Linie durch

den Tropfenschwerpunkt, so ergibt sich ein parabolisches Geschwindigkeitsprofil mit

Scheitelpunkt als absolutem Maximum wcim Tropfenzentrum und dem Minimum wΓ

an der Phasengrenze. Wie bei der L¨

osung von Hadamard-Rybczinski gilt wΓ=−wc.

Das Wirbelzentrum liegt jedoch oberhalb der ¨

Aquatorebene des Tropfens und nicht

innerhalb der ¨

Aquatorebene, wie es bei der L¨

osung von Hadamard-Rybczinski der

Fall ist (siehe Abbildung 7.3). F¨

ur die Tropfen mit λ= 20 und λ= 45 ergibt sich

wc= 0,24 wTund wc= 0,14 wT. Der Betrag der Phasengrenzengeschwindigkeit am

Tropfen¨

aquator ist damit etwa ein Viertel bzw. ein Siebentel der Aufstiegsgeschwin-

digkeit. Im Intervall 20 ≤λ≤45 sinkt die Aufstiegsgeschwindigkeit wTum 0,3 cm/s.

Dies ist einerseits eine signifikante ¨

Anderung, andererseits entspricht dies nur einer

relativen Absenkung von wTum 0,3%. In einem doppelt logarithmischen Diagramm

f¨

ur CDin Abh¨

angigkeit von Re (wie Abbildung 7.2) ist eine solche ¨

Anderung prak-

tisch nicht sichtbar, sodaß Haljasmaas Absch¨

atzung in erster N¨

aherung zutreffend

ist. Aus Abbildung 7.4 ist ersichtlich, daß ab λ≥100 die Aufsiegsgeschwindig-

keit konstant bleibt. Unterhalb λ= 20 steigt die Geschwindigkeit der Tropfen mit

σ= 20.1 mN/m zun¨

achst stark an, was mit einer starken Abflachung des Tropfens ein-

hergeht (Abbildung 7.5). Um den Einfluß von λauf die Aufstiegsgeschwindigkeit wT

ohne Formver¨

anderung studieren zu k¨

onnen, wurde f¨

ur eine weitere Serie von Trop-

fen mit d=5,7 mm und σ= 75 mN/m der Aufstieg in Wasser berechnet. Durch die

7. Simulation von Mais¨

oltropfen

uberh¨

ohte Oberfl¨

achenspannung sind diese Tropfen nahezu kugelf¨

ormig und ¨

Ande-

rungen der Aufstiegsgeschwindigkeiten werden im Wesentlichen durch ¨

Anderungen

von λbestimmt.

0,01 0,1 1 10 100

= 75 mN/m

= 20,1 mN/m

[cm /s]

Abbildung 7.4.: Aufstiegsgeschwindigkeiten f¨

ur Tropfen in Abh¨

angigkeit vom Visko-

sit¨

atenverh¨

altnis λbei d=5,7 mm und ρd= 917 kg/m3. Quadrate: Oberfl¨

achenspannung

σ= 20,1 mN/m. Kreise: σ= 75 mN/m.

0,01 0,1 1 10 100

0,5

0,6

0,7

0,8

0,9

1,0

= 75 mN /m

= 20.1 m N/m

/ d

Abbildung 7.5.: Durchmesserverh¨

altnis f¨

ur Tropfen in Abh¨

angigkeit vom Viskosit¨

aten-

verh¨

altnis λbei d=5,7 mm und ρd= 917 kg/m3. Quadrate: Oberfl¨

achenspannung σ=

20,1 mN/m. Kreise: σ= 75 mN/m.

7.4. Trajektorien

W¨

ahrend des freien Aufstiegs fluider Partikel k¨

onnen im Experiment periodische seit-

liche Oszillationen auftreten, so daß sich eine helikale (d.h. schraubenf¨

ormige) oder

zickzack-f¨

ormige Bewegungsbahn ausbildet. Um zu pr¨

ufen, welche Klassen von Bewe-

gungsbahnen in der Simulation auftreten, wird der Tropfenaufstieg zun¨

achst f¨

ur zwei

exemplarische Tropfengr¨

oßen ohne die Verwendung von Symmetrieebenen bei einem

Wandabstand dW=8dberechnet: Tropfen mit d=5,7 mm steigen im Experiment ge-

radlinig auf, w¨

ahrend Tropfen mit d=13,2 mm horizontal oszillieren. Durch periodi-

sche Randbedingungen an den Seitenw¨

anden wird in der Simulation eine ungehinderte

seitliche Bewegung zugelassen. Bei großen Tropfen stellen sich nach wenigen Sekun-

den freien Aufstiegs periodische seitliche Bewegungen ein, die eher zickzack-f¨

ormig

als helikal sind (siehe Abbildung 7.6). Bei den kleinen Tropfen sind in der Simula-

tion nicht-periodische seitliche Bewegungen festzustellen. Die Trajektorien verlaufen

innerhalb der Ebenen y=x oder y=-x. Weil diese Diagonalen im kubischen Gitter

ausgezeichnete Raumrichtungen darstellen, k¨

onnten numerische Effekte - sogenannte

”parasit¨

are Str¨

omungen“, vgl. Koebe (2004) - die Ursache f¨

ur die unregelm¨

aßige Ho-

rizontalbewegung sein. Um dieses Ph¨

anomen zu unterdr¨

ucken, werden kleine Tropfen

unter Verwendung zweier vertikaler Symmetrieebenen simuliert und so ein geradlini-

ger Aufstieg erzwungen.

Abbildung 7.7 zeigt die zeitliche Entwicklung der Aufstiegsgeschwindigkeit f¨

ur Trop-

fen aus Rechnungen mit zwei Symmetrieebenen. Der kleine Tropfen erreicht nach

einer Beschleunigungsphase eine konstante Endgeschwindigkeit wT. Es pr¨

agt sich ein

station¨

arer Nachlauf aus. Der große Tropfen hingegen wird nach Erreichen einer Ma-

ximalgeschwindigkeit wmax wieder abgebremst. Diese Abbremsung geht einher mit

dem Auftreten von Instabilit¨

aten im Nachlauf, die schließlich zu Wirbelabl¨

osungen

f¨

uhren.

Wird durch Verwendung von weniger als zwei Symmetrieebenen eine seitliche Bewe-

gung des Tropfenschwerpunkts zugelassen, so gehen die Wirbelabl¨

osungen mit einer

zickzack-f¨

ormigen Aufstiegsbewegung einher. Um die Rechenzeit gering zu halten,

wird in diesem Fall nur eine Symmetrieebene verwendet. Dadurch wird die Bewegung

des Tropfenschwerpunkts auf nur eine Ebene festgelegt. Dies schr¨

ankt die Physikalit¨

der Simulation nicht ein, da zickzack-f¨

ormige Trajektorien nur zweidimensional sind.

7. Simulation von Mais¨

oltropfen

-3,5 -3,0 -2, 5 -2,0 -1,5

0,0

0,5

1,0

1,5

2,0

2,5

y [cm]

x [cm ]

Abbildung 7.6.: Trajektorie eines Mais¨

oltropfens (d=13,2 mm) in einer Simulation ohne

Symmetrieebenen (Draufsicht).

0 2 4 6 8 10

d = 13.2 mm

w [c m/s]

t [s ]

d = 5,7 mm

Abbildung 7.7.: Aufstiegsgeschwindigkeit in Abh¨

angigkeit von der Zeit f¨

ur mit zwei Sym-

metrieebenen simulierte Mais¨

oltropfen (d=5,7 mm und d=13,2 mm).

7.4. Trajektorien

Insgesamt wird die Trajektorie eines Tropfens nach folgendem Verfahren eingestuft:

Bildet sich in einer Simulation mit zwei Symmetrieebenen ein stabiler Nachlauf, gilt

der Tropfen als geradlinig aufsteigend. Andernfalls wird die Rechnung mit nur einer

Symmetrieebene wiederholt. Die zeitliche Entwicklung des Geschwindigkeitsfeldes f¨

einen Tropfen mit d=13,2 mm ist im Anhang A.4 dargestellt. Zun¨

achst beschleunigt

der Tropfen geradlinig. W¨

ahrend dieser Phase wird der Nachlauf des Tropfens im-

mer l¨

anger, bis er instabil wird. Schließlich l¨

osen sich Wirbel ab, wechselweise an der

linken und der rechten Tropfenseite. Dies geht einher mit der horizontalen Oszillati-

on. F¨

ur Luftblasen in Wasser wird dieser Vorgang ebenfalls in direkten numerischen

Simulationen beobachtet (Koebe, 2004). Ein detailliertes Modell zur Entstehung der

abwechselnden Wirbelabl¨

osung an den Blasenseiten beschreiben Fan und Tsuchiya

(1990).

In den Abbildungen 7.8 und 7.9 sind die Aufstiegsgeschwindigkeit und die horizon-

tale Position des Schwerpunkts f¨

ur zwei verschieden große Tropfen dargestellt. Beide

Tropfen erreichen bei t=4 sdie gleiche Maximalgeschwindigkeit wmax. Ab t=5 sbe-

ginnt eine periodische horizontale Oszillation. Die Aufstiegsgeschwindigkeit schwankt

dann um einen Mittelwert, der f¨

ur den gr¨

oßeren Tropfen deutlich niedriger ist als f¨

den kleineren Tropfen. F¨

ur oszillierende Tropfen wird eine Endgeschwindigkeit wT

durch diesen Mittelwert definiert. Aus Abbildung 7.9 ist ersichtlich, daß der gr¨

oßere

Tropfen st¨

arker oszilliert. Im Intervall 9 s < t < 10 sverringert sich die seitliche Bewe-

gung des gr¨

oßeren Tropfens kurz. Gleichzeitig steigt die Vertikalgeschwindigkeit. Die

Aufstiegsgeschwindigkeit sinkt demnach mit steigender Amplitude der Oszillation.

7. Simulation von Mais¨

oltropfen

0 2 4 6 8 10 12 14 16

d=13,2

w [c m/s]

t [s ]

d=11,3

Abbildung 7.8.: Aufstiegsgeschwindigkeit in Abh¨

angigkeit von der Zeit f¨

ur mit einer Sym-

metrieebene simulierte Mais¨

oltropfen (d=11,3 mm und d=13,2 mm).

0 2 4 6 8 10 12 14 16

-1,6

-1,4

-1,2

-1,0

-0,8

-0,6

-0,4

-0,2

0,0

0,2

d=11,3 mm

d=13,2 mm

y [c m]

t [s ]

Abbildung 7.9.: Seitliche Verschiebung des Tropfenschwerpunkts in Abh¨

angigkeit von der

Zeit f¨

ur mit einer Symmetrieebene simulierte Mais¨

oltropfen (d=11,3 mm und d=13,2 mm).

7.5. Aufstiegsgeschwindigkeiten und Tropfenformen

In Simulationen mit zwei Symmetrieebenen entwickelt sich f¨

ur Mais¨

oltropfen bis zu

einem Durchmesser von d=9 mm ein stabiler Nachlauf. Demnach werden sie als gerad-

linig aufsteigende Tropfen eingestuft. Bei gr¨

oßeren Tropfen wird der Nachlauf instabil.

Ihre Aufstiegsbewegung wird mit einer Symmetrieebene berechnet, wobei sich eine

zickzack-f¨

ormige Trajektorie einstellt. Die jeweiligen Aufstiegsgeschwindigkeiten wT,

wmax und wosz sind in Abbildung 7.10 wiedergegeben. Im Bereich 9 mm < d < 13 mm

nimmt die Amplitude der seitlichen Oszillation stark zu, was sinkende wosz zur Folge

hat. Ab d=13 mm ist wosz nahezu unabh¨

angig vom Tropfendurchmesser. F¨

ur gerad-

linig aufsteigende Tropfen ist wTin der Simulation etwas kleiner als im Experiment,

f¨

ur oszillierende Tropfen ist wosz ungef¨

ahr gleich groß. Oszillationen treten im Expe-

riment ab d=11mm auf. Die Durchmesserverh¨

altnisse der Tropfen sind in Abbildung

7.11 dargestellt. Je h¨

oher die Aufstiegsgeschwindigkeit ist, desto flacher werden die

Tropfen.

2 4 6 8 10 12 14 16 18 20

Exp.

Sim. lin.

Sim. max .

Sim. osz.

[cm/ s]

d [m m]

Abbildung 7.10.: Aufstiegsgeschwindigkeit f¨

ur Mais¨

oltropfen in Abh¨

angigkeit vom Trop-

fendurchmesser. Kreise: Experiment. Dreiecke: Simulation. Geschlossene Dreiecke: wTf¨

geradlinigen Aufstieg. Offene Dreiecke: wosz. Halboffene Dreiecke: wmax.

7. Simulation von Mais¨

oltropfen

2 4 6 8 10 12 14 16 18 20

0,0

0,1

0,2

0,3

0,4

0,5

0,6

0,7

0,8

0,9

1,0

Exp.

Simula tion:

lin.

max.

osz.

/ d

d [m m]

Abbildung 7.11.: Durchmesserverh¨

altnisse f¨

ur Mais¨

oltropfen in Abh¨

angigkeit vom Trop-

fendurchmesser. Kreise: Experiment. Dreiecke: Simulation. Geschlossene Dreiecke: Gerad-

liniger Aufstieg. Offene Dreiecke: Mittelwert w¨

ahrend seitlicher Oszillation. Halboffene

Dreiecke: Durchmesserverh¨

altnis zum Zeitpunkt des Auftretens der Maximalgeschwindig-

keit wmax .

7.6. Einfluß der tr¨

agen Masse fluider Partikel auf das

Aufstiegsverhalten

Oftmals werden aufsteigende oder sinkende fluide Partikel durch die dimensionslosen

Kennzahlen λ, Re, E¨

o und Mo gekennzeichnet. Die Dichtedifferenz ∆ρist als Faktor

in E¨

o und Mo enthalten. Im Allgemeinen werden E¨

o und Mo auf der Basis von |∆ρ|

berechnet, so daß sich positive Kennzahlen ergeben. Haljasmaa weist darauf hin, daß

hierbei die Information verloren geht, ob es sich um ein aufsteigendes oder sinkendes

Partikel handelt. Bei gleichem |∆ρ|unterscheiden sich aufsteigendes und sinkendes

Partikel durch ihre tr¨

age Masse. Bei gleichf¨

ormiger Bewegung (wie dem geradlini-

gen Aufstieg oder Sinken mit konstanter Geschwindigkeit wT) tritt die tr¨

age Masse

nicht in Erscheinung. W¨

ahrend einer beschleunigten Bewegung jedoch wirkt die tr¨

age

Masse den beschleunigenden Kr¨

aften entgegen. Haljasmaa vergleicht zickzack-f¨

ormi-

ge fluide Partikel mit mechanischen harmonischen Oszillatoren. Sind diese reibung-

behaftet, so sinkt die Amplitude erzwungener Schwingungen mit wachsender tr¨

ager

7.6. Einfluß der tr¨

agen Masse fluider Partikel auf das Aufstiegsverhalten

Masse. Frei aufsteigende Tropfen sind keine harmonischen Oszillatoren, weil es keine

zu irgendeiner Auslenkung proportionale R¨

uckstellkraft gibt. Dennoch ist die An-

nahme plausibel, daß die Amplitude einer horizontalen Oszillation mit wachsender

tr¨

ager Masse des Partikels kleiner wird. Um diese Annahme zu ¨

uberpr¨

ufen, wird ein

Modellfluid untersucht, das bis auf die Dichte gleiche Stoffeigenschaften wie Mais¨

hat. Als Dichte wird ρd= 1097 kg/m3gew¨

ahlt, sodaß die Dichtedifferenz des Modell-

fluids gegen¨

uber Wasser den gleichen Betrag hat wie bei Mais¨

ol (|∆ρ|= 81 kg/m3).

Ein Tropfen des Modellfluids hat wegen seiner gr¨

oßeren Dichte eine um 18% gr¨

oßere

tr¨

age Masse als ein gleich großer Mais¨

oltropfen. Die Auftriebskraft beider Tropfen hat

den gleichen Betrag, jedoch unterschiedliches Vorzeichen. Bei gleicher Gr¨

oße oszilliert

ein Tropfen des Modellfluids in der Simulation weniger stark als ein Mais¨

oltropfen

(Abbildung 7.12), einhergehend mit einer h¨

oheren Aufstiegsgeschwindigkeit der Mo-

delltropfen (Abbildung 7.13). W¨

ahrend eines geradlinigen Aufstiegs/Sinkens sind die

Betr¨

age der Geschwindigkeiten f¨

ur Partikel beider Fluide gleich.

0 2 4 6 8 10 12 14 16 18 20

-1,5

-1,0

-0,5

0,0

0,5

model fluid

y [c m]

t [s]

corn oil

Abbildung 7.12.: Horizontale Verschiebung des Tropfenschwerpunkts in Abh¨

angigkeit von

der Zeit. Mais¨

ol und Modellfluid, d=13,2 mm.

7. Simulation von Mais¨

oltropfen

2 4 6 8 10 12 14 16 18 20

Maisöl

Modell

[cm/ s]

d [m m]

Abbildung 7.13.: Simulierte Aufstiegsgeschwindigkeiten f¨

ur Tropfen in Abh¨

angigkeit vom

Durchmesser. Dreiecke: Mais¨

ol. Sterne: Modellfluid. Geschlossene Symbole: Geradliniger

Aufstieg. Offene Symbole: Horizontale Oszillation. Halboffene Symbole: Maximalgeschwin-

digkeit vor Einsetzen der Oszillation.

8. Rizinus¨

Die dynamische Viskosit¨

at von Rizinus¨

ol ist 900 mal h¨

oher als die von Wasser. Aus

Abbildung 8.1 ist ersichtlich, daß in einem Rizinustropfen die inneren Geschwindig-

keiten vernachl¨

assigbar klein sind.

0.00

0.10

0.20

0.30

0.40

Abbildung 8.1.: Geschwindigkeitsfeld eines Rizinus¨

oltropfens (d=4 mm) nach Erreichen

der Endgeschwindigkeit wT. (Bild um 90◦gedreht, Aufstiegsrichtung nach links.)

8.1. Aufstiegsgeschwindigkeiten und Tropfenformen

Im Experiment steigen Rizinustropfen langsamer auf als Mais¨

oltropfen, was im We-

sentlichen auf die h¨

ohere Dichte des Rizinus¨

ols zur¨

uckzuf¨

uhren ist. Dadurch ist die

Auftriebskraft eines Rizinustropfens nur halb so groß wie die eines Mais¨

oltropfens.

Die Aufstiegsgeschwindigkeit simulierter Mais¨

oltropfen liegt im Falle geradlinigen

Aufstiegs etwas unterhalb der von Haljasmaa (2006) gemessenen Aufstiegsgeschwin-

digkeiten. Im Falle oszillierender Trajektorien sind die Geschwindigkeiten f¨

ur Simu-

lation und Experiment ¨

ahnlich. In der Simulation tritt Wirbelabl¨

osung im Nachlauf

schon ab d=10 mm auf, im Experiment hingegen erst ab d=13,5 mm. Im Bereich

13 mm ≤d≤13,5 mm werden im Experiment sowohl geradlinige als auch zickzack-

f¨

ormige Trajektorien beobachtet. Entsprechend gabelt sich der Graph wT(d) in diesem

Bereich in zwei ¨

Aste auf (englisch: Bifurcation). F¨

ur Tropfendurchmesser d > 13 mm

werden nur zickzack-f¨

ormige Trajektorien beobachtet. In der Simulation tritt das

Ph¨

anomen der Bifurcation nicht auf: F¨

ur eine bestimmte Tropfengr¨

oße ist der Nach-

lauf entweder (d < 10 mm) stabil oder instabil (d≥10 mm).

8. Rizinus¨

0,0 0,5 1,0 1,5 2,0

Exp. lin.

Exp. osz.

Sim. lin.

Sim. osz.

[cm/s]

d [m m]

Abbildung 8.2.: Aufstiegsgeschwindigkeiten f¨

ur Rizinus¨

oltropfen in Abh¨

angigkeit vom

Durchmesser d. Schwarze Symbole: Geradliniger Aufstieg. Graue Symbole: Horizontale Os-

zillation. Kreise: Experiment. Quadrate: Simulation.

0 2 4 6 8 10 12 14 16 18 20

0,5

0,6

0,7

0,8

0,9

1,0

/ d

d [m m]

Abbildung 8.3.: Durchmesserverh¨

altnisse f¨

ur Rizinus¨

oltropfen in Abh¨

angigkeit vom

Durchmesser d. Schwarze Symbole: Geradliniger Aufstieg. Graue Symbole: Horizontale Os-

zillation. Kreise: Experiment. Quadrate: Simulation.

9. Widerstandsbeiwerte f¨

ur ¨

Oltropfen

Basierend auf dem ¨

Aquivalenzdurchmesser dkann jedem Tropfen eine Reynoldszahl

und ein Widerstandsbeiwert zugeordnet werden:

Re = ρdd wT

µd

; (9.1)

CD=4

∆ρ

ρc

g d

wT2.(9.2)

Dabei beruht es lediglich auf Konvention, in (9.1) den ¨

Aquivalenzdurchmesser dzu

benutzen. Es k¨

onnte beispielsweise auch eine Reynoldszahl basierend auf dem hori-

zontalen Tropfendurchmesser dHdefiniert werden. Ebenso wird ¨

uber (9.2) f¨

ur nicht-

sph¨

arische Partikel ein Widerstandsbeiwert definiert, weil diese Gleichung unter der

Voraussetzung einer Kugelform aus (2.1) hergeleitet wurde. Ber¨

ucksichtigt man, daß

f¨

ur nicht-sph¨

arische Partikel der Anstr¨

omquerschnitt durch den Horizontaldurchmes-

ser dHfestgelegt wird, so ergibt sich stattdessen

CD=4

∆ρ

ρc

g d

wT2

.(9.3)

Im allgemeinen sind fluide Partikel abgeflacht (d < dH), weswegen (9.3) kleinere Wi-

derstandsbeiwerte liefert als (9.2). Insbesondere f¨

ur zeitlich nicht formstabile fluide

Partikel ist es jedoch durchaus sinnvoll, als charakteristische L¨

ange den ¨

Aquivalenz-

durchmesser dzu verwenden. In diesem Falle wird durch (9.2) ein Widerstandsbei-

wert bestimmt, der nicht identisch mit dem durch (2.1) definierten Widerstandsbei-

wert ist. Unter Verwendung von (9.1) und (9.2) k¨

onnen Graphen f¨

ur die Aufstiegsge-

schwindigkeit von der dimensionsbehafteten Form wT(d) in die dimensionslose Form

CD(Re) ¨

uberf¨

uhrt werden. Diese Darstellungsform erlaubt auch den Vergleich mit

der Standard-Drag-Kurve (2.6) f¨

ur feste Kugeln. Aus Abbildung 9.1 ist ersichtlich,

daß kleinere ¨

Oltropfen im Experiment ¨

ahnliche Widerstandsbeiwerte aufweisen wie

feste Kugeln. Mit wachsendem Tropfendurchmesser w¨

achst jedoch der Abstand zur

Standard-Drag-Kurve, weil die Tropfen immer flacher werden und daher eine gr¨

oßere

9. Widerstandsbeiwerte f¨

ur ¨

Oltropfen

angestr¨

omte Querschnittsfl¨

ache aufweisen als Kugeln. Ab Re=1200 steigt der Wider-

standsbeiwert der ¨

Oltropfen im Experiment in Verbindung mit oszillierenden Trajek-

torien sprunghaft an. Im Experiment setzt die Oszillation bereits unter Re=1000 ein.

Die Standard-Drag-Kurve gilt nur f¨

ur fixierte angestr¨

omte Kugeln. An frei aufstei-

1000

300

Experi mentel l

Maisöl

Rizinus

Nume risch

Maisöl lin.

Maisöl osz.

Rizinus lin.

Rizinus osz.

3000

0,5

Standard-Dra g

Schirmb lasen: C

= 8/3

Abbildung 9.1.: Widerstandsbeiwerte in Abh¨

angigkeit von der Reynoldszahl f¨

ur ¨

Oltrop-

fen. Gelbe Kreise: Mais¨

oltropfen im Experiment. Schwarze Kreise: Rizinus¨

oltropfen im Ex-

periment. Vielecke: Simulationen. Geschlossene Vielecke: Geradlinige Trajektorien. Offene

Vielecke: Oszillierende Trajektorien. Dreiecke: Mais¨

ol. Quadrate: Rizinus¨

ol.

genden Kugeln beobachtet Haljasmaa ebenfalls einen sprunghaften Anstieg f¨

ur CD,

weil auch sie ab einer bestimmten Gr¨

oße bzw. Reynoldszahl nicht geradlinig aufstei-

gen sondern oszillieren. Auch hier ist dies eine Folge periodischer Wirbelabl¨

osung.

Ausgehend von der Vorstellung, daß die tr¨

age Masse fluider Partikel einer Oszillation

entgegenwirkt, vermutete Haljasmaa zun¨

achst, daß f¨

ur fluide Partikel mit gr¨

oßerer

Dichte als Wasser die kritische Reynoldszahl ¨

uber 1200 liegt. In seinen Experimen-

ten hingegen beobachtete er f¨

ur sinkende fluide Partikel ein Einsetzen der Oszillation

ab Re=800. Dieses Ph¨

anomen f¨

uhrt Haljasmaa auf Formoszillationen (Oszillationen

des Durchmesserverh¨

altnisses) zur¨

uck, die besonders bei Tropfen mit hoher Dichte

auftreten. Formoszillationen l¨

osen m¨

oglicherweise seitliche Oszillationen aus.

10. Zusammenfassung des ersten Teils

Unter Verwendung des massiv parallelen VoF-Codes FS3D wird der Aufstieg von

Tropfen in fl¨

ussigen Medien untersucht. Der Code wird validiert, indem der langsame

Aufstieg von Tropfen in einer sehr viskosen Fl¨

ussigkeit berechnet wird. Die hierf¨

ur nu-

merisch ermittelten Druck- und Geschwindigkeitsfelder stimmen mit der analytischen

L¨

osung f¨

ur den Kriechfall gut ¨

uberein. Unterschiede lassen sich darauf zur¨

uckf¨

uhren,

daß die numerische Simulation auf ein begrenztes Volumen (einen Rechteckkanal)

beschr¨

ankt ist, wohingegen die analytische L¨

osung f¨

ur einen unendlich ausgedehnten

Raum gilt.

F¨

ur kleine, geradlinig aufsteigende Tropfen wird die Aufstiegsgeschwindigkeit in Ab-

h¨

angigkeit von der inneren Viskosit¨

at und der Oberfl¨

achenspannung numerisch un-

tersucht. Je h¨

oher die Viskosit¨

at des Tropfens ist, um so kleiner ist die Aufstiegsge-

schwindigkeit. Wenn die Tropfenviskosit¨

at etwa hundertmal h¨

oher ist als die Visko-

sit¨

at der umgebenden Fl¨

ussigkeit, erreicht die Aufstiegsgeschwindigkeit ihr Minimum.

Bei weiterer Erh¨

ohung der Tropfenviskosit¨

at ¨

andert sich die Aufstiegsgeschwindigkeit

nicht mehr signifikant. Die inneren Geschwindigkeiten im Tropfen werden sehr klein

im Vergleich zu den ¨

außeren Geschwindigkeiten. Die tangentiale Geschwindigkeit der

Phasengrenze wird vernachl¨

assigbar klein, sodaß das fluide Partikel ¨

ahnlich wie ein

festes Partikel gleicher Form umstr¨

omt wird. Tropfen mit hoher Oberfl¨

achenspannung

steigen im Allgemeinen schneller auf, weil sie durch hydrodynamische Effekte nicht

so stark abgeflacht werden und daher eine kleinere angestr¨

omte Querschnittsfl¨

ache

aufweisen.

F¨

ur schnell aufsteigende Tropfen existiert keine analytische L¨

osung. Dies gilt beispiels-

weise f¨

ur den Aufstieg von ¨

Oltropfen in Wasser. F¨

ur diesen Fall werden Experimente

zum Vergleich mit Simulationen herangezogen. Kleine ¨

Oltropfen steigen im Experi-

ment geradlinig auf, w¨

ahrend der Schwerpunkt großer Tropfen w¨

ahrend des Aufstiegs

horziontal oszilliert, woraus sich zickzack-f¨

ormige oder helikale Trajektorien ergeben.

In Simulationen konnte gezeigt werden, dass diese Oszillation durch eine periodische

10. Zusammenfassung des ersten Teils

Wirbelabl¨

osung hervorgerufen wird: Große Tropfen beschleunigen geradlinig bis der

Nachlauf instabil wird und die Wirbelabl¨

osung einsetzt. Diese f¨

uhrt zu einer sponta-

nen Verringerung der Aufstiegsgeschwindigkeit. Da eine periodische Wirbelabl¨

osung

auch hinter festen, in Str¨

omungskan¨

alen fixierten Kugeln beobachtet wird, setzt die

periodische Wirbelabl¨

osung ihrerseits offenbar keine Oszillation des umstr¨

omten Par-

tikels voraus. Wenn das Geschwindigkeitsfeld zum Startzeitpunkt durch Initialisation

mit Zufallswerten gest¨

ort ist, setzt die seitliche Oszillation deutlich fr¨

uher ein. In

der Simulation wird bereits bei kleineren Tropfendurchmessern als im Experiment

nicht-geradliniger Aufstieg beobachtet. Kleine, geradlinig aufsteigende ¨

Oltropfen sind

in der Simulation etwas langsamer als im Experiment. F¨

ur große, nicht geradlinig

aufsteigende ¨

Oltropfen ergibt sich eine befriedigende ¨

Ubereinstimmung f¨

ur die Auf-

stiegsgeschwindigkeit.

Die simulierten Tropfen sind deutlich flacher als im Experiment. Weil im Experiment

Leitungswasser benutzt wurde, k¨

onnten Verunreinigungen (”Surfactants“) hier eine

Rolle gespielt haben. Die Anlagerung von surfactants an eine freie Oberfl¨

ache kann

komplexe Oberfl¨

achenph¨

anomene - sogenannte Marangoni-Effekte - hervorrufen. Um

diese zu unterbinden, m¨

ussten die Experimente in hochreinem Wasser durchgef¨

uhrt

werden. Außerdem wurden nat¨

urliche ¨

Ole verwendet, die ihrerseits verschiedenste

Molek¨

ule enthalten. Um eine homogene, chemisch eindeutig definierte Grenzfl¨

ache zu

erhalten, sollten auch die Tropfen aus reinen Stoffen erzeugt werden.

Teil II.

Blasenaufstieg in linearen Scherstr¨

omungen

11. Stand des Wissens

11.1. Blasen in vertikalen Rohren

Steigen fluide Partikel in ruhenden Fl¨

ussigkeiten auf, so wirken auf sie im Wesent-

lichen Auftriebskraft und Widerstandskraft. Bewegt sich die umgebende Fl¨

ussigkeit

(d.h. die kontinuierliche Phase) in einem durch ¨

außere W¨

ande begrenzten Gebiet,

so treten im Allgemeinen Inhomogenit¨

aten in ihrem Geschwindigkeitsfeld auf. In ei-

ner Rohrleitung zum Beispiel ist die Fließgeschwindigkeit in der Rohrachse maxi-

mal, an den W¨

anden ist sie Null. Durchstr¨

omt eine Fl¨

ussigkeit ein vertikales Rohr in

Aufw¨

artsrichtung, so kann man an darin aufsteigenden Blasen eine seitliche Migra-

tion beobachten: Kleine Blasen wandern w¨

ahrend ihres Aufstiegs in Richtung Rohr-

wand, große Blasen zur Rohrmitte (Zun, 1980, 1988). Detaillierte Untersuchungen zur

Blasenmigration im vertikalen Rohr wurden von Prasser et al. (2005) mittels eines

Gittersensors durchgef¨

uhrt: Ein Gitter aus feinen Dr¨

ahten wird senkrecht zur Haupt-

str¨

omungsrichtung in den Rohrquerschnitt gespannt. An den Kreuzungspunkten der

Dr¨

ahte wird simultan die Leitf¨

ahigkeit in der Str¨

omung gemessen. Passiert eine Bla-

se einen Kreuzungspunkt, ist die Leitf¨

ahigkeit nahezu Null. Ist der Kreuzungspunkt

vollst¨

andig mit Wasser umgeben, so wird hohe Leitf¨

ahigkeit gemessen. So kann f¨

jeden Zeitpunkt die Verteilung des Gasgehaltes ¨

uber den Rohrquerschnitt gemessen

werden. Gr¨

oßere Blasen legen gleichzeitig mehrere benachbarte Gitterknoten trocken

und k¨

onnen so individuell registriert und vermessen werden. Auf diese Weise ist es

m¨

oglich, die radiale Verteilung des Gasgehaltes in einem Rohr nach Gr¨

oßenklassen

aufzuschl¨

usseln. In Abbildung 11.1 sind solche Verteilungen f¨

ur einen Rohrradius von

97,5 mm dargestellt. F¨

ur das System Wasser/Luft unter Normalbedingungen erkennt

man eine deutlich erh¨

ohte Konzentration kleiner Blasen in der N¨

ahe der Rohrwand.

Beispielsweise haben Blasen mit einem Durchmesser von 0 −3,5 mm in einer Entfer-

nung von 70 mm von der Rohrachse nur einen Volumenanteil von knapp 0,4%. Bei

95 mm ist ihr Volumenanteil doppelt so hoch. Im Gegensatz hierzu f¨

allt der Volumen-

11. Stand des Wissens

anteil der Blasen mit einem Durchmesser von 6,5−7 mm zum Rand hin stark ab.

Im System Wasser/Dampf bei 6,5 M Pa zeigt nur die Gr¨

oßenklasse bis 0 −3,5 mm

einen deutliches Maximum an der Rohrwand.

Abbildung 11.1.: Radiale Verteilung des Gasgehaltes in einem Rohr mit einem Radius

von 97,5 mm, aufgeschl¨

usselt nach Blasengr¨

oßenklassen (¨

Aquivalenzdurchmesser in mm).

Links: Luftblasen in Wasser unter Normalbedingungen. Rechts: Dampfblasen in siedendem

Wasser bei 6,5 M Pa.

Ist die Fl¨

ussigkeitsstr¨

omung in einem vertikalen Rohr abw¨

arts gerichtet, so kehrt sich

der Effekt um: Die kleinen Blasen wandern zur Rohrmitte und die großen nach au-

ßen. Dieses Ph¨

anomen deutet darauf hin, daß auf Blasen in inhomogenen str¨

omenden

Fl¨

ussigkeiten zus¨

atzlich zu Auftriebs- und Widerstandskraft eine weitere Kraft wirkt,

die Liftkraft genannt wird. Eine genaue Kenntnis dieser Kraft ist erforderlich, wenn

die Blasenbewegung in chemischen Reaktoren korrekt simuliert werden soll - etwa

mit dem Euler-Euler oder dem Euler-Lagrange Verfahren. Hier wird ein Schließungs-

term ben¨

otigt, der die Liftkraft in Abh¨

angigkeit von den Eigenschaften der Blase

und den Geschwindigkeitsfeldern beider Phasen modelliert. So wurde beispielsweise

durch Lucas et al. (2005) gezeigt, daß die Liftkraft einen erheblichen Einfluß auf den

Betriebszustand von Blasens¨

aulen hat. Die Liftkraft kann dazu beitragen, daß der

Str¨

omungszustand in Blasens¨

aulen instation¨

ar wird.

11.1. Blasen in vertikalen Rohren

Abbildung 11.2.: Geschwindigkeitsfeld einer linearen Scherstr¨

omung.

Der einfachste Fall einer inhomogenen Str¨

omung ist die lineare Scherstr¨

omung, wie

sie in Abb. 11.2 dargestellt ist. Hier nimmt die vertikale Geschwindigkeitskomponente

proportional zur Verschiebung in einer horizontalen Richtung zu. Die Proportiona-

lit¨

atskonstante ωhat die Einheit [1/s] und wird ”Scherrate“ genannt.

vl=





−ωx





⇒rot vl=











=: ω(11.1)

Die Rotation eines Geschwindigkeitsfeldes ist eine vektorielle Gr¨

oße, die angibt, mit

welcher Winkelgeschwindigkeit ein Lagrange’sches Partikel im Str¨

omungsfeld rotie-

ren w¨

urde. F¨

ur die in Abbildung 11.2 dargestellte Scherstr¨

omung ist anschaulich klar,

daß ein infinitesimal kleines, festes Partikel im Uhrzeigersinn rotieren w¨

urde, weil es

auf der ”rechten“ Seite st¨

arker angestr¨

omt wird. Also zeigt ωin die Bildebene hinein

(Rechte-Hand-Regel). In einer linearen Scherstr¨

omung ist die Rotation ωr¨

aumlich

und zeitlich konstant und soll ”Scherratenvektor“ genannt werden. Es gilt |ω|=ω.

Bewegt sich ein Partikel mit einer Relativgeschwindigkeit vrin dieser Scherstr¨

omung,

so erf¨

ahrt es eine Liftkraft senkrecht zu vrund zu ω. Nach umfassenden Experimenten

zur Blasenbewegung in einem aufw¨

arts durchstr¨

omten vertikalen Rohr mit quadrati-

schem Querschnitt schl¨

agt Zun (1980) folgenden Term f¨

ur die Liftkraft vor:

FL=−CLVpρlvr×ω.(11.2)

Hier ist Vpdas Partikelvolumen, ρldie Dichte der Fl¨

ussigkeit und CLein dimensions-

loser Koeffizient. Die Relativgeschwindigkeit ist die Differenz der Partikelgeschwin-

digkeit vpund der Str¨

omungsgeschwindigkeit vl, wie sie im ungest¨

orten Fall am Ort

11. Stand des Wissens

des Partikelschwerpunkts vorl¨

age. Gleichung (11.2) kann als Schließungsterm f¨

ur die

Liftkraft verwendet werden, wenn der Wert von CLbekannt ist.

11.2. Theoretische Untersuchungen zu kugelf¨

ormigen Blasen in

linearen Scherstr¨

omungen

F¨

ur feste, kugelf¨

ormige Partikel in linearen Scherstr¨

omungen konnte Auton (1987)

auf analytischem Wege den Liftkoeffizienten bestimmen. Unter der Annahme eines

idealen Fluids und einer ”kleinen“ Scherrate konnte er zeigen, daß CL= 0,5 gilt.

”Klein“ heißt hier, daß die ¨

Anderung von vLentlang des Partikeldurchmessers D

klein gegen¨

uber der Relativgeschwindigkeit ist:

ωD ≪ |vr|mit vr=vp−vl,

oder anders formuliert

Sr ≪1 mit Sr = ωD

|vr|,

wobei Sr eine dimensionslose Scherrate ist. In Mehrphasenstr¨

omungen mit Scheranteil

tritt sie als f¨

unfte dimensionslose Kennzahl zu Re, Mo, Eo und dem Dichteverh¨

alt-

nis hinzu. Unter Verweis auf Autons Resultat f¨

ur feste Kugeln wird in der Litera-

tur CL= 0,5 oftmals auch f¨

ur kugelf¨

ormige Blasen angenommen (Legrende und

Magnaudet, 1998; Tomiyama, 2004). Feste Kugeln weisen an der Phasengrenze eine

Haftbedingung auf. Dies gilt auch f¨

ur Gasblasen, wenn die Oberfl¨

ache vollst¨

andig mit

Surfactants belegt ist. F¨

ur Gasblasen in vollkommen reinem Wasser gilt an der Ober-

fl¨

ache Schlupfbedingung, wenn die Viskosit¨

at des Gases vernachl¨

assigt wird. In der

Herleitung von Auton ist die Fl¨

ussigkeit ideal (µl= 0). In diesem Fall geht die Haft-

bedingung in eine Schlupfbedingung ¨

uber. Die Anwendung von Autons Resultat auf

Kugelblasen setzt voraus, daß die Viskosit¨

at der Fl¨

ussigkeit ”hinreichend“ klein ist,

um die Annahme einer idealen Fl¨

ussigkeit zu rechtfertigen. Ist die Blasenoberfl¨

ache

ungleichm¨

aßig stark mit Surfactants belegt, so variert die Oberfl¨

achenspannung und

es treten sogenannte Marangoni-Effekte auf (Alke, 2007). Derartige Effekte sind hier

nicht ber¨

ucksichtigt.

F¨

ur sph¨

arische Blasen in reibungsbehafteten linearen Scherstr¨

omungen haben Legren-

de und Magnaudet (1998) umfassende theoretische Untersuchungen durchgef¨

uhrt, die

11.3. Das Experiment von Tomiyama

von einer Kugel mit Schlupfbedingung als Modell f¨

ur eine kugelf¨

ormige Blase ausge-

hen. Als Resultat analytischer Ans¨

atze in Verbindung mit numerischen Rechnungen

erhalten sie eine Korrelation f¨

ur den Liftkoeffizienten in Abh¨

angigkeit von Reynolds-

zahl und Scherrate:

CL=s6

π2

2.255

(Re Sr)1/2(1 + 0.2 Re/Sr)3/22

+1

1 + 16/Re

1 + 29/Re2

.(11.3)

Im Falle Re ≫Sr ist der Nenner des ersten Terms proportional zu Re2/Sr, im Falle

Re ≪Sr ist er proportional zu (Re Sr)1/2. Im Grenzfall Re → ∞ verschwindet bei

fester Scherrate der erste Term und es ergibt sich insgesamt CL= 1/2. Im Grenzfall

Re →0 wird CLunendlich groß. Die Abh¨

angigkeit des Liftkoeffizienten von der

Scherrate wird im Bereich Re >5 verschwindend gering.

11.3. Das Experiment von Tomiyama

Die prominentesten Experimente zu der Bewegung von Blasen in linearen Scher-

str¨

omungen wurden von Tomiyama und Mitarbeitern durchgef¨

uhrt (Tomiyama et

al., 1995, 2002; Tomiyama, 2004). Der Versuchsaufbau besteht aus einem hohen qua-

derf¨

ormigen Tank, der mit einer Wasser-Glycerin-L¨

osung gef¨

ullt ist. In dem Tank

verl¨

auft parallel zu einer Wand in senkrechter Richtung ein breites Gummiband, wel-

ches einem F¨

orderband ¨

ahnelt. In dem 30 mm breiten Spalt zwischen Wand und

Band entsteht eine lineare Scherstr¨

omung, sobald das Band in Bewegung gesetzt

wird. Die Geschwindigkeit des Bandes legt die Scherrate im Spalt fest. Durch eine

D¨

use werden Luftblasen in den Spalt injiziert, deren Form und Bewegung mit Hoch-

geschwindigkeitskameras aufgezeichnet wird. Der Tank wird mehrfach mit Wasser

und Glycerin bef¨

ullt, jeweils in unterschiedlichen Mischungsverh¨

altnissen. Daraus er-

geben sich unterschiedliche Viskosit¨

aten bei nur geringen ¨

Anderungen der Dichte und

der Oberfl¨

achenspannung (siehe Tabelle 12.1). Tomiyama et al. (2002) zeichneten die

Aufstiegsbahnen verschieden großer Blasen (d= 3...6 mm) in L¨

osungen mit Visko-

sit¨

aten von 18 bis 89 mPa s auf. Sie verwendeten destilliertes Wasser, um den Einfluß

von Verunreinigungen (surfactants) zu minimieren. Je nach Gr¨

oße migrierten die Bla-

sen w¨

ahrend ihres Aufstiegs im Spalt entweder zur Tankwand oder zum F¨

orderband.

Unter Verwendung der Bewegungsgleichung (12.6) f¨

ur Blasen wurden die jeweiligen

Liftkoeffizienten mittels eines nicht n¨

aher beschriebenen numerischen Verfahrens aus

11. Stand des Wissens

den Trajektorien ermittelt. F¨

ur kleine Blasen wurden positive Liftkoeffizienten er-

mittelt, deren Betr¨

age von der Reynoldszahl abh¨

angig sind. F¨

ur mittlere bis große

Blasen k¨

onnen Vorzeichen und Betrag als Funktion einer modifizierten E¨

otv¨

os-Zahl

E¨

oHbeschrieben werden, die sich auf den horizontalen Durchmesser dHder Blase

bezieht:

E¨

oH=g∆ρ d2

σ.(11.4)

Insgesamt faßt Tomiyama die Ergebnisse des Experiments in folgender Korrelation

zusammen:

CL=(min[f1(Re), f2(E¨

oH)] f¨

ur E¨

oH<4

f2(E¨

oH) f¨

ur E¨

oH≥4,(11.5)

wobei

f1(Re) = 0,288 tanh(0,121 Re) (11.6)

und

f2(E¨

oH) = 0,00105 E¨

H−0,0159 E¨

H−0,00204 E¨

oH+ 0,474 .(11.7)

Diese Korrelation ist rein empirisch und wurde durch Interpolation der Datenpunk-

te f¨

ur Blasen im Bereich Re > 5 erhalten. Demnach wechselt CLbei E¨

oH≈6 das

Vorzeichen. F¨

ur gr¨

oßere Blasen ist CLnegativ, d.h. sie migrieren in entgegengesetzte

Richtung verglichen mit kleineren Blasen. Durch die theoretischen Resultate von Le-

grende und Magnaudet (1998) motiviert, wurde das Experiment f¨

ur kleine Blasen

(Re <5) wiederholt (Tomiyama, 2004). Abh¨

angig von der Viskosit¨

at der Fl¨

ussig-

keit waren die Ergebnisse sehr unterschiedlich. Bei µl= 55 mPa s wurden hin zu

kleinen Reynoldszahlen sehr große positive Liftkoeffizienten (bis zu CL= 10) festge-

stellt, in ¨

Ubereinstimmung mit Gleichung (11.3). Bei kleineren Viskosit¨

aten konnten

jedoch auch negative Liftkoeffizienten (bis zu CL=−2) beobachtet werden, was im

Widerspruch zu Gleichung (11.3) steht.

12. Simulation von Blasen in Glycerinl¨

osungen

Mit dem Code FS3D ist es m¨

oglich, die Bewegung von Blasen in Scherstr¨

omungen

zu berechnen, wobei die Verformung der Phasengrenze Bestandteil der numerischen

L¨

osung ist. Umfassende Simulation mit FS3D zum Blasenaufstieg in ruhenden w¨

assri-

gen Systemen wurden durch Koebe (2004) durchgef¨

uhrt.

12.1. Randbedingungen

F¨

ur verschiedene Wasser-Glycerin L¨

osungen und Blasendurchmesser wird das Expe-

riment von Tomiyama nachgerechnet. Die Stoffwerte der Fluide sind Tabelle 12.1

zu entnehmen, die Blasendurchmesser gehen aus Abb. 12.2 hervor. Testrechnungen

ergeben, daß ein Wandabstand von vier Blasendurchessern in y-Richtung und von

acht Blasendurchmessern in x-Richtung (die Richtung der Migration) ausreicht, um

die Wandeinfl¨

usse auf die Blasenbewegung zu eliminieren. Weil die Randbedingungen

in y-Richtung vollkommen symmetrisch sind, ist eine Spiegelsymmetrie der L¨

osung

bez¨

uglich der x-y-Ebene zu erwarten. Daher kann das Rechengebiet durch Verwen-

dung einer Spiegelebene am Rand um die H¨

alfte verkleinert werden. Eine Spiegel-

symmetrie als Randbedingung bedeutet eine Schlupfbedingung f¨

ur das Geschwindig-

keitsfeld, und einen Kontaktwinkel von 90◦f¨

ur die Phasengrenze, d.h. die Phasen-

grenze muß senkrecht auf die Spiegelebene auftreffen. Um eine lineare Scherstr¨

omung

in der Fl¨

ussigphase im Rechengebiet zu gew¨

ahrleisten, sind die beiden W¨

ande, die

sich in x-Richtung gegen¨

uberliegen, mit einer Haftbedingung ausgestattet. Die Wand

bei x= 0 ist fest, die Wand bei x= 8 dbewegt sich mit einer festen Geschwin-

digkeit U0=−ω8d. Demnach sind an diesen W¨

anden die Geschwindigkeiten mit

v= (0,0,0) bzw. v= (0,0, U0) fest vorgegeben (Dirichlet’sche Randbedingung).

Als Standardaufl¨

osung werden 32 Zellen pro Blasendurchmesser gew¨

ahlt. Am oberen

Rand des Rechengebietes (bei z=0) wird v=−ωx vorgegeben, d.h. Fl¨

ussigphase

wird von oben mit einer bereits voll ausgepr¨

agten Scherstr¨

omung eingeleitet. Die

12. Simulation von Blasen in Glycerinl¨

osungen

gegeneinander bewegten W¨

ande erhalten dieses Geschwindigkeitsprofil im gesamten

Rechengebiet aufrecht.

= xω

bewegte

Wand

bewegte Wand

Seitansicht Draufsicht

Ausstroembedingung

Einstroemung

feste

Wand

feste Wand

8 d

Schlupfbedingung

Symmerie− −ebene

2 d

Abbildung 12.1.: Geometrie des Rechengebietes. Randbedingungen f¨

ur die Simulation

von Blasen mit FS3D

Tabelle 12.1.: Stoffdaten f¨

ur Wasser-Glycerin L¨

osungen aus (Tomiyama et al., 2002) (µl=

18...53 mPa s). Zum Vergleich sind auch die Stoffdaten von Wasser aufgef¨

uhrt (µl= 1 mPa s)

µl[mPa s] 1 18 19 53

ρl[kg/m3] 1000 1150 1154 1180

σ[mN/m] 72 63 61 63

log10(Mo) -10,6 -5,5 -5,3 -3,6

12.2. Blasenformen und Aufstiegsgeschwindigkeiten

Das experimentelle Ergebnis von Tomiyama et al. (2002) l¨

aßt vermuten, daß zwischen

Blasenform und Liftkraft eine Wechselbeziehung besteht, denn die E¨

otv¨

os-Zahl geht

maßgeblich in die empirische Korrelation ein. Die Blasenform wiederum ist eng mit

12.2. Blasenformen und Aufstiegsgeschwindigkeiten

der E¨

otv¨

os-Zahl verkn¨

upft: Blasen mit kleiner E¨

otv¨

os-Zahl sind eher kugelf¨

ormig, so-

chle mit großer E¨

otv¨

os-Zahl st¨

arker deformiert (abgeflacht). Dieser Zusammenhang

wird ersichtlich aus der Bedeutung der E¨

otv¨

os-Zahl als Quotient aus Auftriebskraft

durch die Kraft der Oberfl¨

achenspannung (vgl. Kap. 5.2). Dominiert die Oberfl¨

achen-

spannung, so ist die Blase eher kugelf¨

ormig. Dominiert die Auftriebskraft, so wird die

Blase st¨

arker durch hydrodynamische Effekte deformiert und insbesondere abgeflacht.

Darauf wird in Kapitel 13.2.2 noch n¨

aher eingegangen.

2 3 4 5 6

19 mPa s EXP

19 mPa s SIM

53 mPa s EXP

53 mPa s SIM

[mm]

d [m m ]

Abbildung 12.2.: Horizontaler Blasendurchmesser dHin Abh¨

angigkeit vom ¨

Aquivalenz-

durchmesser d. Schwarze Symbole: µl= 53 mPa s. Graue Symbole: µl= 19 mPa s. Offene

Symbole: Blasenformen im Experiment. Geschlossene Symbole: Simulierte Blasenformen.

Durchgezogene Linie: Kugelform (dH=d).

In Abbildung 12.2 werden die experimentell beobachteten mit den simulierten Bla-

senformen verglichen. Die ¨

Ubereinstimmung ist sehr gut: Die simulierten Horizon-

taldurchmesser weichen maximal um 2% von den beobachteten ab. Die durchgezo-

gene Linie repr¨

asentiert kugelf¨

ormige Blasen (dH=d). Je weiter Datenpunkte ¨

uber

dieser Linie liegen, desto st¨

arker sind die zugeh¨

origen Blasen abgeflacht. Blasen in

der weniger viskosen L¨

osung sind breiter (und flacher) als in der viskoseren L¨

osung.

Große Blasen sind flacher als kleine. Die simulierten Horizontaldurchmesser sind un-

12. Simulation von Blasen in Glycerinl¨

osungen

abh¨

angig von der Scherrate; insbesonders sind die Horizontaldurchmesser in x- wie

in y-Richtung stets gleich, trotz Anisotropie des Geschwindigkeitsfeldes der Fl¨

ussig-

phase. Die Abflachung der Blasen h¨

angt eng mit ihrer Aufstiegsgeschwindigkeit zu-

sammen: Schnelle Blasen sind flacher als langsame. Daher liegt es nahe, auch die

simulierten Aufstiegsgeschwindigkeiten mit experimentellen Werten zu vergleichen.

Tomiyama et al. (2002) schlagen drei verschiedene semi-empirische Korrelationen f¨

den Widerstandsbeiwert als Funktion der Reynoldszahl vor - je nach Verunreinigung

der Fl¨

ussigphase:

F¨

ur reine Systeme

CD= max min 16

Re(1 + 0.15Re0.687),48

Re,8

E¨

o + 4; (12.1)

f¨

ur leicht kontaminierte Systeme

CD= max min 24

Re(1 + 0.15Re0.687),72

Re,8

E¨

o + 4; (12.2)

f¨

ur vollst¨

andig kontaminierte Systeme

CD= max 24

Re(1 + 0.15Re0.687),8

E¨

o + 4.(12.3)

Im Grenzfall Re →0 geht die Korrelation f¨

ur reine Systeme in die Beziehung CD=

16/Re ¨

uber, die sich aus der L¨

osung von Hadamard (1911) und Rybczynski (1911)

f¨

ur Blasen ergibt. Die Korrelationen f¨

ur kontaminierte Systeme gehen in CD= 24/Re

uber, was der Stoke’schen L¨

osung f¨

ur feste Kugeln entspricht. In der konkreten An-

wendung der Korrelation (12.1) auf die simulierten Blasen wird der Dragkoeffizient

entweder durch das erste Argument der Minimumsfunktion (f¨

ur kleine Blasen) oder

durch die E¨

otv¨

oszahl (f¨

ur große Blasen) festgelegt. Mittels eines iterativen Algorith-

mus k¨

onnen die Beziehungen von der Form CD(Re,E¨

o) in die Form wT(d)¨

uberf¨

uhrt

werden. In Abbildung 12.3 werden die simulierten Aufstiegsgeschwindigkeiten mit

den Resultaten der Korrelationen verglichen. F¨

ur µl= 53 mPa s sind auch die Auf-

stiegsgeschwindigkeiten in kontaminiertem Wasser eingezeichnet. Im entprechenden

Wertebereich vereinfachen sich die Korrelationen (12.2) und (12.3) zu

CD=24

Re (1 + 0,15Re0.687).(12.4)

12.3. Aufstiegsbahnen

Dies ist eine empirische Korrelation f¨

ur den Widerstandsbeiwert starrer Kugeln. Ins-

besonders f¨

ur kleine Blasen stimmen die simulierten Aufstiegsgeschwindigkeiten bes-

ser mit der Korrelation f¨

ur reines Wasser ¨

uberein. Dies ist zu erwarten, denn eine

Belegung der Phasengrenze mit Surfactants ist in der Simulation nicht modelliert.

12.3. Aufstiegsbahnen

Die Aufstiegsgeschwindigkeiten aus Abbildung 12.3 sind in ruhender Fl¨

ussigkeit be-

rechnet worden, d.h. bei ω= 0. Steigt eine Blase im Abstand x von der festen Wand im

Rechengebiet auf, so erf¨

ahrt sie eine Gegenstr¨

omung der Fl¨

ussigphase mit dem Betrag

ωx (siehe Abb. 11.2). Ihre absolute Aufstiegsgeschwindigkeit vpim ortsfesten Bezugs-

system h¨

angt also von x ab. Eine ortsunabh¨

angige Aufstiegsgeschwindigkeit kann in

Scherstr¨

omungen nur ¨

uber die Relativgeschwindigkeit definiert werden. Hierzu bildet

man die Differenz der absoluten Aufstiegsgeschwindigkeit vpund der Fl¨

ussigkeitsge-

schwindigkeit vl, wie sie im ungest¨

orten Fall (d.h. ohne Anwesenheit der Blase) am

Ort des Blasenschwerpunkts vorl¨

age. F¨

ur die lineare Scherstr¨

omung ergibt sich

vr=vp−vl=vp−(0,0,−ωx).(12.5)

Im Folgenden sollen u,vund wjeweils die Komponenten der Relativgeschwindigkeit

vrbezeichnen. Wegen der Symmetrieebene bei den Randbedingungen (vgl. Abbildung

12.1) wird eine Blasenbewegung in y-Richtung in den Simulationen nicht zugelassen,

d.h. vist stets Null. Die Zeitabh¨

angigkeit von uund wist in den Abbildungen 12.4 und

12.5 f¨

ur zwei exemplarische Blasen dargestellt. Zum Zeitpunkt t= 0 starten jeweils

kugelf¨

ormige Blasen in einer linearen Scherstr¨

omung. Bereits ab etwa 0,2 s stellen

sich konstante Geschwindigkeiten in vertikaler Richtung ein (Abb. 12.4). Die relative

Endgeschwindigkeit win vertikaler Richtung ist jeweils identisch mit den absoluten

Endgeschwindigkeiten wTgleich großer Blasen in ruhender Fl¨

ussigkeit. Auch in hori-

zontaler Richtung stellt sich f¨

ur alle Blasen eine konstante Endgeschwindigkeit uein,

mit Ausnahme der großen Blase bei µl= 53 mPa s, bei der uged¨

ampft oszilliert (Abb.

12.5). In diesem Fall l¨

aßt sich durch Mittelung der Extrema ein konstanter Endwert

f¨

ur uextrapolieren. In ¨

Ubereinstimmung mit dem Experiment kommen Oszillatio-

nen der Horizontalgeschwindigkeit unur bei großen Blasen und geringen Viskosit¨

aten

µlvor. In den ersten hundertstel Sekunden ist umaximal (besonders deutlich bei

12. Simulation von Blasen in Glycerinl¨

osungen

2 3 4 5 6

19 mPa s

53 mPa s

[cm/s]

d [mm]

c o n te r m i n a te d

cl ea n

Abbildung 12.3.: Aufstiegsgeschwindigkeiten wTin Abh¨

angigkeit vom ¨

Aquivalenzdurch-

messer d. Datenpunkte: Simulierte Werte. Durchgezogene Linien: Korrelation (12.1) aus-

gewertet f¨

ur µl= 19 mPa s (schwarz) und f¨

ur µl= 53 mPa s (grau). Gestrichelte Linie:

Korrelation (12.2) f¨

ur kontaminiertes Wasser, ausgewertet f¨

ur µl= 53 mPa s.

d=2,8 mm). Dies geht einher mit einer starken Abflachung der Blasen unmittelbar

nach dem Start als Kugelform. Wenn die Blasen ihre endg¨

ultige stabile Form er-

reicht haben, stimmen ihre Horizontaldurchmesser dHmit denen aus Abbildung 12.2

uberein. Es wird also kein Einfluß der Scherrate auf dHbeobachtet. In allen untersuch-

ten F¨

allen ist die dimensionslose Scherrate Sr klein (Sr <0,26). Bei kleinen Blasen

stellt sich eine positive Horizontalgeschwindigkeit uein, d.h. sie migrieren in ein Ge-

biet mit st¨

arkerer Gegenstr¨

omung, wodurch ihre absolute Aufstiegsgeschwindigkeit

sinkt. Die großen Blasen migrieren seitlich in Richtung Gegenstr¨

omung (u < 0) mit

steigender absoluter Aufstiegsgeschwindigkeit. Eine Bewegung mit konstanten Re-

lativgeschwindigkeiten einer linearen Scherstr¨

omung bedeutet im ortsfesten Bezugs-

system eine gleichf¨

ormige Horizontalbewegung und eine gleichm¨

aßig beschleunigte

Vertikalbewegung. Daraus ergibt sich (wie beim schiefen Wurf) eine parabelf¨

ormige

Bewegungsbahn.

12.3. Aufstiegsbahnen

0,0 0,1 0, 2 0,3

d = 5,6 mm

d = 3,2 mm

d = 2,8 mm

w [c m/s]

t [s ]

Abbildung 12.4.: Vertikale Komponente der Relativgeschwindigkeit simulierter Blasen in

Abh¨

angigkeit von der Zeit. Graue Kurven: µl= 19 mPa s, ω= 8.3s−1. Schwarze Kurven:

µl= 19 mPa s, ω= 6.2s−1

0,0 0,1 0, 2 0,3 0,4 0,5

d = 2,8 mm

d = 3,2 mm

d = 5,6 mm

u [cm /s]

t [s ]

Abbildung 12.5.: Horitzontale Komponente der Geschwindigkeit simulierter Blasen.

12. Simulation von Blasen in Glycerinl¨

osungen

12.4. Die Bestimmung des Liftkoeffzienten CLaus

Blasentrajektorien

Betrachtet man die Blase als Massepunkt, so lautet ihre Bewegungsgleichung

(ρg+ 0,5ρl)dvp

dt=−3CDρl

4d|vr|vr−CLρlvr×rotvl+ (ρg−ρl)g.(12.6)

Diese Gleichung ist eine Gleichung f¨

ur Kraftdichten. Wenn man sie mit dem Bla-

senvolumen VBmultipliziert, so erh¨

alt man die Gleichung f¨

ur Kr¨

afte. Die Terme von

links nach rechts sind: Tr¨

agheitskraft, Widerstandskraft, Liftkraft und Auftriebskraft.

F¨

ur die tr¨

age Masse der Blase wird VB(ρg+ 0,5ρl) angenommen. Es wird also da-

von ausgegangen, daß die Blase zus¨

atzlich zur Masse des Gases in ihrer unmittelba-

ren Umgebung (beispielsweise in der N¨

ahe der Phasengrenze oder im Nachlauf) eine

Wassermenge mitf¨

uhrt, die der H¨

alfte ihres Volumens entspricht. Die Masse 0,5VBρg

wird virtuelle Masse genannt. Gegen¨

uber der virtuellen Masse ist die Masse des Gases

vern¨

achl¨

assigbar klein. Der Faktor 0,5 ist lediglich ein Sch¨

atzwert. In der folgenden

Rechnung wird die virtuelle Masse eliminiert, weshalb das Problem ihrer Absch¨

atzung

hier ohne Bedeutung ist.

Dividiert man die Bewegungsgleichung auf beiden Seiten durch die Fl¨

ussigkeitsdichte

ρlund vernachl¨

assigt man die Gasdichte (ρg≪ρl), so vereinfacht sie sich zu

0,5dvp

dt=−D|vr|vr−CLvr×ω−gmit D=3CD

4d.(12.7)

Zun¨

achst soll die linke Seite berechnet werden. F¨

ur die Lagrange’sche Ableitung einer

skalaren Gr¨

oße Φ entlang einer Blasentrajektorie gilt:

d Φ

dt=∂tΦ + (vp·∇)Φ .

Anstelle von Φ wird vp=vl+vreingesetzt. Weil vlund vrnicht explizit von der

Zeit abh¨

angen und nur vlortsabh¨

angig ist, bleibt

dvp

dt= (vp·∇)vl.

Nach Einsetzen von vrund vlergibt sich

dvp

dt=−





ωu







.(12.8)

12.4. Die Bestimmung des Liftkoeffzienten CLaus Blasentrajektorien

Dies bedeutet, daß die Blase sich in horizontaler Richtung mit gleichbleibender Ge-

schwindigkeit bewegt. In der vertikalen Richtung erf¨

ahrt sie jedoch eine gleichm¨

aßige

Beschleunigung ω u . Durch Einsetzen von Gleichung (12.8) und (11.1) in (12.7) ergibt

sich unter der Annahme einer konstanten Relativgeschwindigkeit vr= (u, 0, w)

−0,5ω











=−D|vr|









−CLω





−w







+











.(12.9)

Hier ist g=|g|positiv, weil die Auftriebskraft in Richtung +z wirkt. Gl. (12.9)

kann als Gleichungssystem mit einer Gleichung f¨

ur jede Raumrichtung interpretiert

werden. Die zweite Gleichung f¨

ur die y-Richtung ist trivial (0 = 0). F¨

ur die anderen

beiden Raumrichtungen ergeben sich die Gleichungen

0 = −D|vr|u+CLωw f¨

ur die x-Richtung, (12.10)

und

−0,5ω u =−D|vr|w−CLωu +gf¨

ur die z-Richtung . (12.11)

Dieses Gleichungssystem enth¨

alt die beiden Unbekannten CLund D, wobei letztere

den Dragkoeffizienten CDenth¨

alt. Durch Aufl¨

osen von (12.10) nach Dund Einsetzen

in (12.11) kann Deliminiert werden:

−0,5ωu =−CLωw2

u−CLωu +g . (12.12)

Aufl¨

osen nach CLliefert:

CL=u

u2+w2g

ω+ 0,5u.(12.13)

F¨

ur die Simulationen gilt stets 0,5u≪g/ω und u2≪w2. Somit kann man (12.13)

vereinfachen zu:

CL=g u

ω w2.(12.14)

12. Simulation von Blasen in Glycerinl¨

osungen

12.5. Die Abh¨

angigkeit des Liftkoeffizienten von der Scherrate

F¨

ur Blasen mit den ¨

Aquivalenzdurchmessern d=3,56 mm und d=4,85 mm werden

die Aufstiegsbahnen in der L¨

osung mit µl= 18 mPa s (vgl. Tabelle 12.1) f¨

ur ver-

schiedene Scherraten simuliert und jeweils der Liftkoeffizient nach Gleichung (12.14)

ermittelt. Aus Abbildung 12.8 ist ersichtlich, daß der Liftkoeffizient von der Scherrate

unabh¨

angig ist. Dies steht nicht im Widerspruch zum Ergebnis von Magnaudet (vgl.

Gleichung 11.3), weil die hier untersuchten Blasen im Bereich mittlerer Reynolds-

zahlen liegen (Re >5). Die Vorzeichen der simulierten Liftkoeffizienten stimmen mit

den experimentellen Beobachtungen von Tomiyama et al. (2002) ¨

uberein. F¨

ur den

Blasendurchmesser d=4,85 mm f¨

allt jedoch auf, daß der Betrag des simulierten Lift-

koeffizienten etwas gr¨

oßer ist als im Experiment. Im Folgenden wird stets auf das

Experiment von Tomiyama et al. (2002) Bezug genommen. Wo nicht ausdr¨

ucklich

erw¨

ahnt, wird f¨

ur Simulationen die Standardaufl¨

osung von 32 Zellen pro ¨

Aquivalenz-

durchmesser dverwendet.

1 2 3 4 5 6 7 8 9

-0,3

-0,2

-0,1

0,0

0,1

0,2

0,3

d = 3,56 mm

d = 4,85 mm

[1/s]

Abbildung 12.6.: Liftkoeffizient f¨

ur Blasen (d=3,56 und d=4,85) in Abh¨

angigkeit von der

Scherrate (µl= 18 mPa s). Offene Symbole: Experiment. Geschlossene Symbole: Simulation.

F¨

ur ω= 8s−1liegt kein experimenteller Wert vor.

12.6. CLin Abh¨

angigkeit vom Blasendurchmesser

12.6. CLin Abh¨

angigkeit vom Blasendurchmesser

F¨

ur Blasen verschiedener Gr¨

oßen (d=2,8-5,6 mm) wird die Bewegung in Glycerinl¨

osun-

gen bei jeweils konstant gehaltener Scherrate simuliert und die Liftkoeffizienten in

Abh¨

angigkeit von der modifizierten E¨

otv¨

os-Zahl E¨

oHdargestellt. Fast alle simulier-

ten Werte liegen unterhalb der Kurve, welche die experimentelle Korrelation (11.5)

repr¨

asentiert. Zudem zeigt sich in der Simulation - anders als im Experiment - eine

deutliche Abh¨

angigkeit der simulierten Werte von der Viskosit¨

at der Fl¨

ussigphase ab.

Es k¨

onnen verschiedene Ursachen f¨

ur die Unterschiede zwischen Simulation und Ex-

periment verantwortlich sein:

a) Die Aufl¨

osung des Gitters ist bei 32 Zellen/dm¨

oglicherweise noch nicht hoch ge-

nug, um alle relevanten Str¨

omungsph¨

anomene abzubilden.

b) Im L¨

oser treten numerische Fehler auf (beispielsweise in Form unphysikalischer

Geschwindigkeitsfelder, sogenannter ”parasit¨

arer Str¨

omungen“).

c) Nicht alle im Experiment auftretenden Ph¨

anomene sind in dem mathematischen

Modell ber¨

ucksichtigt, das dem numerischen L¨

oser zugrunde liegt. Beispielsweise ist

der Einfluß von Surfactants nicht ber¨

ucksichtigt. Zwar wurde im Experiment destil-

liertes Wasser zur Herstellung der L¨

osungen verwendet, aber dennoch k¨

onnten geringe

Verunreinigungen eine Rolle spielen.

0 2 4 6 8 10

-0,6

-0,4

-0,2

0,0

0,2

0,4

Tom iyama

19 m Pa s

53 m Pa s

Eö

Abbildung 12.7.: Liftkoeffizient f¨

ur Blasen in Abh¨

angigkeit von E¨

oHbei einer Aufl¨

osung

von 32 Zellen/d. Dreiecke: µl= 19 mPa s, bei ω= 6,2s−1. Kreise: µl= 53 mPa s, bei

ω= 8,3s−1.

12. Simulation von Blasen in Glycerinl¨

osungen

12.7. Abh¨

angigkeit der simulierten Liftkoeffizienten von der

Aufl¨

osung

Die Rechnungen aus dem vorangegangenen Abschnitt werden mit verschiedenen Auf-

l¨

osungen (16 bis 64 Zellen pro Blasendurchmesser) wiederholt. Wegen des hohen Re-

chenaufwandes konnte nur ein Fall mit der Aufl¨

osung 64 Zellen/dgerechnet werden.

Die Betr¨

age f¨

ur die Liftkoeffizienten werden f¨

ur große Blasen kleiner und kommen

somit den Werten der Korrelation von Tomiyama n¨

aher. Es ist jedoch nicht zu er-

warten, daß bei weiterer Erh¨

ohung der Aufl¨

osung noch signifikante ¨

Anderungen der

Resultate auftreten, denn bereits bei einer Aufl¨

osung von 32 Zellen/dist die Breite

der Zellen kleiner als die Kolmogorov-L¨

angenskala λK. Diese liefert eine Absch¨

atzung

f¨

ur die Gr¨

oße der kleinsten m¨

oglichen Wirbel in der Fl¨

ussigphase:

λK=νl3

ǫ

(12.15)

mit der kinematischen Viskosit¨

at νl=µl/ρlund der auf die Masse bezogenen Dissi-

pationsrate ǫ. F¨

ur aufsteigende Partikel kann diese abgesch¨

atzt werden als ǫ=g w.

2 4 6 8 10

-0,8

-0,6

-0,4

-0,2

0,0

0,2

0,4

Tomi yama

19 mPa s

16 Zellen / d

32 Zellen / d

64 Zellen / d

53 mPa s

16 Zellen / d

32 Zellen / d

Eö

Abbildung 12.8.: Liftkoeffizient f¨

ur Blasen in Abh¨

angigkeit von E¨

oHbei verschiedenen

Aufl¨

osungen .

12.7. Abh¨

angigkeit der simulierten Liftkoeffizienten von der Aufl¨

osung

Aus den Geschwindigkeitsfeldern ist ersichtlich, daß sich bei µl= 53 mPa s im Nach-

lauf der Blasen keine Wirbel bilden (Abb. 12.9). F¨

ur µl= 19 mPa s sehen die Ge-

schwindigkeitsfelder qualitativ ¨

ahnlich aus. Die Blasen sind lediglich etwas flacher.

Wegen der einfachen Struktur dieser Geschwindigkeitsfelder kann davon ausgegangen

werden, daß alle relevanten str¨

omungsmechanischen Ph¨

anomene bei den verwendeten

Gitteraufl¨

osungen voll erfaßt sind. Es nicht zu erwarten, daß bei h¨

oheren Aufl¨

osun-

gen komplexere Strukturen in den Geschwindigkeitsfeldern - etwa kleine Wirbel -

erkennbar w¨

urden.

Abbildung 12.9.: Geschwindigkeitsfelder um Blasen bei µl= 53 mPa s und ω= 8,3s−1.

Fette Pfeile: Bewegungsrichtung der Blase. Links: Kleine Blase; d=3,2 mm. Rechts: Große

Blase; d=5,6 mm. Die L¨

angeneinheit der Skalierung ist [cm]. Man beachte, daß der Maßstab

der beiden Abbildungen unterschiedlich ist. Die Geschwindigkeit der Gegenstr¨

omung steigt

von links nach rechts.

12. Simulation von Blasen in Glycerinl¨

osungen

Parasti¨

are Str¨

omungen treten bei Zweiphasen-Simulationen an der Phasengrenze auf

(Koebe, 2004). Nach Augenschein sind die Geschwindigkeitsfelder in Abbildungen

12.9 und 12.10 auch an der Phasengrenze physikalisch plausibel. Erfahrungsgem¨

aß

spielen parasit¨

are Str¨

omungen eine Rolle, wenn die gesamte kinetische Energie nahe

Null ist (Scardovelli und Zaleski, 1999). Dies trifft auf die betrachteten Situationen

nicht zu.

Auch in DNS-Simulationen mit anderen Codes sind die berechneten Liftkoeffizienten

systematisch kleiner als nach der Korrelation von Tomiyama (Sousa et al., 2006;

Rusche, 2002; Kuipers et al., 2006). Das spricht f¨

ur die Annahme, daß im Experiment

außer der reinen Zweiphasenstr¨

omung weitere Effekte - wie etwa Surfactants - eine

Rolle spielen.

1.00 1.10 1.20 1.30 1.40 1.50 1.60

1.00

1.10

1.20

1.30

1.40

vof−function

[−]

vmin = 0.551725

vmax = 23.3500

velocity

[cm/s]

t = 0.100001 s

c = 3081

2.0 2.5 3.0 3.5

1.6

1.8

2.0

2.2

2.4

2.6

2.8

Abbildung 12.10.: Geschwindigkeitsfelder um Blasen bei µl= 19 mPa s und ω= 6,2s−1.

Links: Kleine Blase; d=2,8 mm. Rechts: Große Blase; d=5,6 mm. Die L¨

angeneinheit der

Skalierung ist [cm]. Man beachte, daß der Maßstab der beiden Abbildungen unterschiedlich

ist. Die Geschwindigkeit der Gegenstr¨

omung steigt von links nach rechts.

12.8. Variation der Oberfl¨

achenspannung

12.8. Variation der Oberfl¨

achenspannung

Gem¨

aß der Korrelation (11.5) ist der Liftkoeffizient eine Funktion der modifizierten

E¨

otv¨

os-Zahl E¨

oH, mit

E¨

oH=g∆ρ dH2

σ.

Im Experiment ebenso wie in den vorangegangenen Abschnitten wurde E¨

oHdurch die

Variation des Blasendurchmessers ver¨

andert. Prinzipiell kann E¨

oHjedoch auch ¨

uber

die Auftriebskraft oder die Oberfl¨

achenspannung variiert werden. Steigt beispielswei-

se die Oberfl¨

achenspannung σ, so sinkt E¨

oHaus zwei Gr¨

unden: Erstens tritt σim

Nenner von E¨

oHauf, zweitens sinkt der horizontale Blasendurchmesser dH, weil die

Oberfl¨

achenspannung einer Abflachung der Blasen entgegenwirkt. Im Folgenden wer-

den Simulationsergebnisse f¨

ur Serien von Modellsystemen vorgestellt, in denen jeweils

nur die Oberfl¨

achenspannung unterschiedlich ist (Tabelle 12.2). Indem σvariiert wird,

ver¨

andert sich auch die Mortonzahl des Systems, weil im Nenner der Mortonzahl σ

in der dritten Potenz auftritt:

Mo = g∆ρ µ4

ρ2

lσ3

Hebt man etwa in System A die Oberfl¨

achenspannung von 10 auf 75 mN/m, so

sinkt die Mortonzahl um knapp drei Gr¨

oßenordnungen. Derartige Serien von Stoff-

systemen lassen sich in der Natur nicht finden, denn im Allgemeinen unterscheiden

sich unterschiedliche Stoffsysteme meist in mehreren Parametern. Zudem sind die

Oberfl¨

achenspannungen im Modellsystem B unrealistisch hoch. In einer Simulation

hingegen k¨

onnen auch ”unrealistische“ Stoffsysteme untersucht werden. Hier wird bei-

spielsweise der Einfluß der Oberfl¨

achenspannung auf den Liftkoeffizienten untersucht.

Tabelle 12.2.: Serien von Modellsystemen zur Variation der Oberfl¨

achenspannung. Se-

rie A:σ= 10, 15, 20, 30, 40, 50, 60, 75 mN/m. Serie B: σ= 100, 200, 300, 400, 600,

800 mN/m.

Serie d [mm] µl[mPa s] ρl[kg/m3]σ[mN/m] log(Mo)

A 3,2 19 1,154 10...75 -2,9... -5,6

B10 100 1,000 100...800 -3,0... -5,7

12. Simulation von Blasen in Glycerinl¨

osungen

In Abbildung 12.11 ist die Abh¨

angigkeit des Liftkoeffizienten von E¨

oHdargestellt,

wobei E¨

oHauf zweierlei Weise variiert wurde: Durch Simulation von Blasen mit

unterschiedlichem Durchmesser d(bei gleichbleibender Oberfl¨

achenspannung) oder

durch Simulation von Blasen mit verschiedenen Oberfl¨

achenspannungen (bei gleich-

bleibendem Blasendurchmesser). Alle Datenreihen k¨

onnen gut approximiert werden

durch die lineare Funktion

CL=−0,1 E¨

oH+ 0,5.(12.16)

Die Steigung dieser Funktion ist ¨

ahnlich wie die Steigung der Tomiyama-Korrelation

im Bereich 4 <E¨

oH<7. Die lineare Approximation wechselt ihr Vorzeichen bei

E¨

oH= 5, die Korrelation jedoch erst bei ungef¨

ahr E¨

oH≈6. Im Grenzfall unendlich

großer Oberfl¨

achenspannung werden Blasen kugelf¨

ormig und E¨

oHgeht gegen Null.

F¨

ur E¨

oH= 0 liefert Gleichung (12.16) den Wert CL= 0,5 . Dies entspricht genau dem

Liftkoeffizienten, den Auton f¨

ur feste Kugel in linearen Scherstr¨

omungen berechnet

hat.

0 2 4 6 8 10 12 14 16

-1,0

-0,8

-0,6

-0,4

-0,2

0,0

0,2

0,4

Variation of d

Tomiyama

19 mPas

Variation of

19 mPas

100 mPas

y=0,5-0,1x

Abbildung 12.11.: Liftkoeffizient in Abh¨

angigkeit von E¨

oH. Schwarze Kurve: Tomiyama-

Korrelation. Geschlossene Dreiecke: µl= 19 mPa s; σ= 61 mN/m; d=2,8...5,5 mm.

Offene Dreiecke: µl= 19 mPa s; σ= 10...75 mN/m; d=3,2 mm (”Serie A“).

Sterne: µl= 100 mPa s; σ= 100...800 mN/m; d=10 mm (”Serie B“).

12.8. Variation der Oberfl¨

achenspannung

In den Serien A und B werden jeweils fast drei Gr¨

oßenordnungen von Mo ¨

uberstrichen,

wohingegen bei der dritten Datenreihe in Abbildung 12.11 (Variation von dbei µl=

19 mPa s) Mo konstant ist (log Mo=-5,3). Dennoch lassen sich alle Datenreihen durch

eine Funktion von E¨

oHgut approximieren (Abb. 12.11). Offenbar hat die Mortonzahl

des Systems kaum einen Einfluß auf die Liftkraft.

In den Abbildungen 12.12 und 12.13 ist der Liftkoeffizient in Abh¨

angigkeit von der auf

den ¨

Aquivalenzdurchmesser bezogenen E¨

otv¨

oszahl E¨

o beziehungsweise in Abh¨

angig-

keit vom Durchmesserverh¨

altnis dV/dHabgetragen. In beiden Darstellungen haben

die Datenreihen f¨

ur unterschiedliche Systeme jeweils einen anderen Verlauf. Daraus ist

ersichtlich, daß sich CLnicht als Funktion von E¨

o oder dV/dHbeschreiben l¨

aßt. Damit

verbleibt allein E¨

oHals entscheidende Kennzahl, die den Betrag des Liftkoeffizienten

bestimmt. F¨

ur stark abgeflachte Blasen (kleines dV/dH) ist der Liftkoeffizient negativ,

f¨

ur ann¨

ahernd kugelf¨

ormige Blasen (dV/dH≈1) ist er positiv. Sowohl in der Simu-

lation (Abb. 12.13) als auch im Experiment (Abb. 12.14) k¨

onnen bei dV/dH≈0,6

sowohl positive als auch negative Liftkoeffizienten auftreten - je nach Stoffsystem.

0 2 4 6 8 10 12

-1,2

-1,0

-0,8

-0,6

-0,4

-0,2

0,0

0,2

0,4

Various d

19 mPa s

Various

19 mPa s

100 mPa s

y = 0,5 - 0,125 x

Abbildung 12.12.: Liftkoeffizient in Abh¨

angigkeit von E¨

o (auf den ¨

Aquivalenzdurchmesser

dbezogene E¨

otv¨

oszahl). Stoffsysteme wie in Abb. 12.11 und Tabelle 12.2.

12. Simulation von Blasen in Glycerinl¨

osungen

0,4 0,5 0,6 0,7 0,8

-1,2

-1,0

-0,8

-0,6

-0,4

-0,2

0,0

0,2

0,4

Variation of d

19 mPa s

Variation of

19 mPa s

100 mPa s

Abbildung 12.13.: Liftkoeffizient in Abh¨

angigkeit vom Durchmesserverh¨

altnis dV/dH.

Stoffsysteme wie in Abb. 12.11 und Tabelle 12.2.

0,4 0,6 0, 8 1,0

-0,3

-0,2

-0,1

0,0

0,1

0,2

0,3

53 mPa s

19 mPa s

Abbildung 12.14.: Liftkoeffizient in Abh¨

angigkeit vom Durchmesserverh¨

altnis dV/dHf¨

Blasen aus dem Experiment von Tomiyama et al. (2002).

12.9. Variation der Gravitation

Da die E¨

otv¨

oszahl das Verh¨

altnis von Auftriebskraft zur Oberfl¨

achenspannungskraft

angibt, kann E¨

oHauch ¨

uber die Gravitation variiert werden. F¨

ur einen Blase mit

fest vorgegebenem Durchmesser (d=5,5 mm) wird die Trajektorie f¨

ur den Aufstieg

bei unterschiedlicher Gravitation berechnet (g= 9,81; 8; 6; 4; 2 m2/s). Je geringer

die Gravitation, desto langsamer steigt die Blase auf, und die Deformation durch

hydrodynamische Effekte (Abflachung) nimmt ab. Als Fl¨

ussigphase wurde das System

µl= 19 mPa s aus Tabelle 12.1 gew¨

ahlt, mit der Scherrate ω= 6,2s−1. Aus

Abbildung 12.15 ist ersichtlich, daß der Liftkoeffizient in ¨

ahnlicher Weise von E¨

abh¨

angt, wie der Durchmesser der Blasen bei konstanter Gravitation (g= 9,81 m2/s)

variiert wird. In beiden F¨

allen lassen sich alle Werte mit der linaren Funktion (12.16)

approximieren.

0 2 4 6 8 10 12

-0,6

-0,4

-0,2

0,0

0,2

0,4

Tomi yama

d varied

g varied

y =0,5 -0,1x

Eö

Abbildung 12.15.: Liftkoeffizient in Abh¨

angigkeit von E¨

oH. Schwarze Dreiecke: Variation

des Blasendurchmessers (vgl. Abb. 12.7). Graue Kreise: Variation der Gravitation.

12. Simulation von Blasen in Glycerinl¨

osungen

12.10. Luftblasen in Wasser

Von besonderem Interesse ist die Gr¨

oße der Liftkraft auf Blasen in Wasser. Hierf¨

gibt es zwar eine Anzahl von Experimenten, die eine Absch¨

atzung des kritischen Bla-

sendurchmessers erm¨

oglichen, bei dem der Liftkoeffizient sein Vorzeichen wechselt

(beispielsweise Beobachtungen der Blasenverteilung in vertikalen Rohren, vgl. Ab-

schnitt 11.1). Es sind jedoch bisher noch nicht die Betr¨

age des Liftkoeffizienten f¨

Blasen in Wasser experimentell ermittelt worden. Dies liegt unter anderem daran, daß

es wegen der geringen dynamischen Viskosit¨

at schwer ist, in Wasser eine stabile linea-

re Scherstr¨

omung zu erzeugen. Aufgrund der Kelvin-Helmholz-Instabilit¨

at kann bei

h¨

oheren Reynoldszahlen eine lineare Scherstr¨

omung zusammenbrechen und es k¨

onnen

Turbulenzen entstehen. In der Simulation mit FS3D bleibt die Scherstr¨

omung auch f¨

Wasser stabil. Anders als im Experiment kann in der Simulation eine v¨

ollig ungest¨

orte

Scherstr¨

omung als Anfangsbedingung eingestellt werden. St¨

orungen wie thermische

Konvektion oder Ersch¨

utterungen k¨

onnen hier ausgeschlossen werden. Zur Stabilisie-

rung des simulierten Geschwindigkeitsprofils tr¨

agt die Einstr¨

om-Randbedingung mit

einer linearen Scherstr¨

omung am oberen Rand des Rechengebietes bei. Es werden

Trajektorien f¨

ur Blasen (d=1...6 mm) in Wasser berechnet. Die relative Aufstiegsge-

0,0 0,1 0, 2 0,3 0,4 0,5 0, 6 0,7 0,8

-1,0

-0,8

-0,6

-0,4

-0,2

0,0

x = -1,20 t - 0,07

x [cm]

t [s]

x = -1,14 t + 0,08

Abbildung 12.16.: Graue Kurve: Horizontale Koordinate des Schwerpunkts einer Blase

(d=6mm) in Wasser in Abh¨

angigkeit von der Zeit. ω= 10s−1.

12.10. Luftblasen in Wasser

schwindigkeit dieser Blasen liegt zwischen 20 und 24 cm/s, und ist damit vergleichbar

mit experimentell beobachteten Aufstiegsgeschwindigkeiten von Blasen in ruhendem

Wasser (Simulationen von Blasen in ruhenden Fl¨

ussigkeiten und Vergleiche mit Expe-

rimenten hat Koebe (2004) detailliert durchgef¨

uhrt). Anders als in den viskosen Gly-

cerinl¨

osungen oszillieren gr¨

oßere Blasen in Wasser so stark, daß eine zickzack-f¨

ormige

Trajektorie mit periodischer Richtungsumkehr der Horizontalbewegung entsteht. In

diesem Fall stellt sich zwar keine konstante Horizontalgeschwindigkeit ein, jedoch

nach wenigen horizontalen Schwingungen ist die mittlere Horizontalgeschwindigkeit

(gemittelt ¨

uber jeweils eine Schwingungsperiode) konstant. Im Weg-Zeit-Diagramm

k¨

onnen die relativen Maxima oder Minima des Graphen x(t) durch Geraden verbun-

den werden. Die Steigung dieser Geraden ergibt jeweils einen Wert f¨

ur die mittlere Ho-

rizontalgeschwindigkeit. Der Mittelwert der beiden Geradensteigungen wird als Wert

f¨

ur die Horizontalgeschwindigkeit verwendet. Tomiyama wendet die aus dem Expe-

0 2 4 6 8 10

-0,6

-0,4

-0,2

0,0

0,2

0,4

Tom iyama

53 m Pa s

19 m Pa s

1 m Pa s

Eö

Abbildung 12.17.: Liftkoeffizient von Blasen in Abh¨

angigkeit von der modifizierten

E¨

otv¨

os-Zahl E¨

oH. Schwarze Kurve: Korrelation von Tomiyama. Sterne: Simulierte Blasen in

Wasser (µl= 1 mPa s), bei einer Aufl¨

osung von 16 Zellen/d. Schwarze Dreiecke und graue

Kreise: Simulierte Blasen in Glycerinl¨

osungen µl= 53 bzw. 19 mPa s bei einer Aufl¨

osung

von 32 Zellen/d.

12. Simulation von Blasen in Glycerinl¨

osungen

riment mit Glycerinl¨

osungen gewonnene Korrelation auch auf Luftblasen in Wasser

an. Um den Liftkoeffizienten f¨

ur eine bestimmte Blasengr¨

oße dzu ermitteln, ist die

Kenntnis des dazugeh¨

origen Horizontaldurchmessers dHerforderlich. Tomiyama ver-

wendet die Korrelation von Wellek et al. f¨

ur Blasen in vollst¨

andig kontaminiertem

Wasser, um den Horizontaldurchmesser der Blasen abzusch¨

atzen:

1 + 0,163 Eo0,757 .(12.17)

Unter der zus¨

atzlichen Annahme elliptischer Blasenformen (d3=dVd2

H) kann mittels

der Wellek-Korrelation zu jedem ¨

Aquivalenzdurchmesser dder zugeh¨

orige Horizon-

taldurchmesser dHin kontaminiertem Wasser dHermittelt werden. Aus dHund den

Stoffdaten ergibt sich E¨

oH, und ¨

uber Tomiyamas Korrelation (11.5) wird schließlich

der zugeh¨

orige Liftkoeffizient ermittelt. So kann der Liftkoeffizient als Funktion vom

Aquivalenzdurchmesser dargestellt werden werden. Diese Funktion ist in Abbildung

12.19 eingezeichnet. Ihr Nulldurchgang liegt bei d=5,8 mm und stimmt recht gut mit

dem kritischen Blasendurchmesser aus experimentellen Beobachtungen ¨

uberein. Nach

der gleichen Vorgehensweise hat Prasser et al. (2005) den kritischen Durchmesser f¨

Dampfblasen in siedendem Wasser unter einem Druck von 6,5 MPa errechnet. Die

Siedetemperatur bei diesem Druck betr¨

agt 280 K. Wegen der geringen Oberfl¨

achen-

spannung (σ=) ergibt sich aus der Tomiyama-Korrelation ein kritischer Blasendurch-

messer von 3,5 mm, was mit dem Experiment ¨

ubereinstimmt (vgl. Abb. 11.1). Den-

noch ist die Verkn¨

upfung von Wellek-Korrelation (12.17) und Tomiyama-Korrelation

(11.5) nicht konsistent, denn Tomiyama verwendete in seinem Experiment destillier-

tes Wasser, um den Einfluß von Surfactants auszuschließen. Die Wellek-Korrelation

hingegen gilt nur f¨

ur vollst¨

andig kontaminiertes Wasser. Blasen in reinem Wasser

hingegen k¨

onnen deutlich flacher sein. So schlagen Fan und Tsuchiya (1990) f¨

ur Bla-

sen in reinen Systemen eine Korrelation vor, die auf der sogenannten Tadaki-Zahl

Ta = ReMo0,23 beruht:

={0,81 + 0,2 tanh[1,8(0,4−log Ta)]}3f¨

ur 0,3≤Ta ≤20 .(12.18)

Diese Korrelation gilt f¨

ur ellipsoide Blasen. Unterhalb Ta = 0,3 gelte dV/dH= 1.0

(sph¨

arische Blasen), oberhalb Ta = 20 sei dV/dH= 0.24 (Kappenblasen). Aus Ab-

bildung 12.18 geht hervor, daß die Korrelation (12.18) die simulierten Blasenformen

12.10. Luftblasen in Wasser

weitaus besser trifft als diejenige von Wellek. In der Simulation oszilliert der Horizon-

taldurchmesser gr¨

oßerer Blasen in Verbindung mit einer zickzack-f¨

ormigen Aufstiegs-

bahn. In solchen F¨

allen werden Mittelwerte f¨

ur dVund dHverwendet.

01234567

Simula tion

Fan

W el lek

Kugel

[ m m ]

d [m m]

Abbildung 12.18.: Horizontaler Durchmesser dHals Funktion des ¨

Aquivalenzdurchmes-

sers d. Gestrichelte Gerade: Kugelform (dH=d). Schwarze Kurve: Korrellation (12.17) von

Wellek. Graue Kurve: Korrelation (Fan und Tsuchiya (1990)). Schwarze Quadrate: Simu-

lierte Blasen.

12. Simulation von Blasen in Glycerinl¨

osungen

0 2 4 6 8

-0,6

-0,4

-0,2

0,0

0,2

0,4

1 m Pa s

Tom iyama

d [mm]

Abbildung 12.19.: Liftkoeffizient von Blasen in Wasser in Abh¨

angigkeit vom ¨

Aquivalenz-

durchmesser d. Schwarze Kurve: Korrelation von Tomiyama. Sterne: Simulierte Blasen, bei

einer Aufl¨

osung von 16 Zellen/d.

13. Untersuchung der Geschwindigkeits- und

Druckfelder an der Phasengrenze

Alle Kr¨

afte, die zwischen Blase und Fl¨

ussigkeit wirken, werden ¨

uber die Phasengrenze

ausgetauscht. Um das Zustandekommen der Liftkraft (und insbesondere ihres Vor-

zeichenwechsels) zu verstehen, liegt es daher nahe, die Druck- und Geschwindigkeits-

felder an der Phasengrenze zu untersuchen. Es bietet sich an, mit dem zweidimen-

sionalen Fall zu beginnen, weil hier die Felder besonders einfach visualisiert werden

k¨

onnen. Allerdings gibt es zu zweidimensionalen Blasen kein ¨

Aquivalent in der Natur.

Eine kreisf¨

ormige 2D-Blase entspr¨

ache beispielsweise einem unendlich ausgedehntem

Zylinder.

13.1. Zweidimensionale Blasen

Erste numerische Simulationen zu Blasen in Scherstr¨

omungen wurden von Tomiyama

und Magnaudet durchgef¨

uhrt. Hierbei beschr¨

ankte man sich wegen der begrenzten

Rechenleistung auf die Untersuchung zweidimensionaler Blasen. Es konnte gezeigt

werden, daß bereits im zweidimensionalen Fall kleine Blasen in Richtung st¨

arkerer

Gegenstr¨

omung migrieren, wogegen große Blasen in die entgegengesetzte Richtung

wandern. Dieses Ph¨

anomen wird auch in Simulationen mit FS3D beobachtet. In Ab-

bildung 13.1 sind die Geschwindigkeitsfelder um 2D-Blasen in einer Glycerin-¨

ahn-

lichen Modellfl¨

ussigkeit mit ρl= 1,2kg/m3,µl= 100 mPa s und σ= 64 mN/m

dargestellt. Nur f¨

ur die große Blase (d=6 mm) stellt sich nach einer beschleunigten

Startphase eine konstante Relativgeschwindigkeit ein. Die horizontale Migrationsge-

schwindigkeit der kleinen Blase (d=3 mm) erreicht ihr Maximum nach der Startphase,

wenn die maximale Relativgeschwindigkeit in vertikaler Richtung erreicht ist. Danach

f¨

allt die Migrationsgeschwindigkeit jedoch monoton. Diese Abbremsung kann durch

Vergr¨

oßerung des horizontalen Wandabstandes dWverringert werden. Allerdings ist

auch f¨

ur dW= 16dzu keinem Zeitintervall die Horizontalgeschwindigkeit der kleinen

13. Untersuchung der Geschwindigkeits- und Druckfelder an der Phasengrenze

0.20

0.30

0.40

0.50

0.60

0.70

0.80

0.5

1.0

1.5

2.0

Abbildung 13.1.: Geschwindigkeitsfelder um 2D-Blasen bei ω= 10s−1. Links: Kleine

Blase; d=3 mm. Rechts: Große Blase; d=6 mm. Man beachte, daß der Maßstab der beiden

Abbildungen unterschiedlich ist. Der Geschwindigkeit der Gegenstr¨

omung steigt von links

nach rechts.

Blase konstant, weshalb die Bestimmung des Liftkoeffizienten nach Gleichung (12.14)

in diesem Fall unm¨

oglich ist. Zum st¨

arkeren Wandeinfluß in 2D-Simulationen k¨

onnte

folgender Zusammenhang beitragen: Wenn sich eine zweidimensionale Blase in eine

Raumrichtung bewegt, so kann die verdr¨

angte Fl¨

ussigkeit nur in eine andere Rich-

tung ausweichen. Bei einer dreidimensonalen Blase kann verdr¨

angte Fl¨

ussigkeit stets

in zwei Raumrichtungen ausweichen.

13.1.1. Zirkulation um zweidimensionale Blasen

In der kleinen Blase aus Abbildung 13.1 (links) ist der rechte Wirbel ausgepr¨

agter, weil

er von der vorbeistreichenden Fl¨

ussigphase st¨

arker angetrieben wird. In der großen

Blase hingegen ist aufgrund der Deformation der linke Wirbel ausgepr¨

agter. Beide

Blasen migrieren in Richtung des st¨

arker ausgepr¨

agten Wirbels. Dieses Ph¨

anomen

beschreiben auch Ervin und Tryggvason (1997) und vermuten, daß die Richtung der

Migration durch die Zirkulation des Geschwindigkeitsfeldes um die Phasengrenze be-

stimmt wird. Die Zirkulation Zum die Blase ist das Linienintegral

Z=IΓ

vds(13.1)

13.1. Zweidimensionale Blasen

entlang der Phasengrenze Γ gegen den Uhrzeigersinn um die Blase. Dabei ist dsein

infinitesimales Linienelement der Phasengrenze. Anstatt des Skalarprodukts vdskann

auch auch das Kreuzprodukt v×nintegriert werden, mit nals ¨

außerer Normale der

Phasengrenze.

Z=eyIΓ

v×nds . (13.2)

Hierbei ist ds die L¨

ange eines infinitesimalen Linienelementes der Phasengrenze, und

eyein Einheitsvektor senkrecht zur Aufstiegsrichtung und senkrecht zur Richtung der

seitlichen Migration (also senkrecht auf der simulierten 2D-Ebene). Das Linienintegral

(13.2) kann unter Verwendung von n=−∇f/ |∇f|in ein Fl¨

achenintegral ¨

uber das

Rechengebiet G umgewandelt werden:

Z=ZG

ey[v×(−∇f)] dA . (13.3)

In dieser Form wird die Zirkulation um die in Abbildung 13.1 dargestellten Blasen

numerisch berechnet. Es ergibt sich f¨

ur die kleine Blase Z = -0.80 cm2/s und f¨

ur die

große Blase Z=2,37. Teilt man die Zirkulation um eine 2D-Blase durch ihre Fl¨

ache,

so ergibt sich ein Wert f¨

ur die durchschnittliche Rotation der Gasphase:

ωg=Z

A.(13.4)

Auf diese Weise erh¨

alt man f¨

ur die kleine Blase ωg= 8.4s−1und f¨

ur die große Blase

ωg=−11.5s−1. Die Rotation der ungest¨

orten Scherstr¨

omung ist hier ω= 10s−1. In

der großen Blase hat die durchschnittliche Rotation also ein anderes Vorzeichen als

in der Fl¨

ussigphase und in der kleinen Blase. Diese Fallbeispiele st¨

utzen demnach die

Hypothese von Magnaudet, wonach die Richtung der Blasenmigration vom Vorzeichen

der Gaszirkulation abh¨

angt.

13.1.2. Dynamischer Druck um zweidimensionale Blasen

Der Druck pin der Fl¨

ussigphase l¨

aßt sich zerlegen in einen hydrostatischen Druck

ps=p0−ρlgz und einen hydrodynamischen Druckanteil pd:

pd:= p−ps=p−p0+ρlgz .

13. Untersuchung der Geschwindigkeits- und Druckfelder an der Phasengrenze

Abbildung 13.2.: Kleine Blase, d=3 mm (vgl. Abb. 13.1) Links: Feld des dynamischen

Drucks in der Umgebung der Blase. Rechts: Dynamischer Druck nahe der Phasengrenze als

Funktion der H¨

ohe h¨

uber der Blasenbasis. Graue Quadrate: Dynamischer Druck auf die

linke Blasenseite. Graue Rauten: Dynamischer Druck auf die rechte Blasenseite. Schwarze

Kreise: Differenz der beiden Dr¨

ucke.

Abbildung 13.3.: Große Blase, d=6 mm (vgl. Abb. 13.1) Links: Feld des dynamischen

Drucks in der Umgebung der Blase. Rechts: Dynamischer Druck nahe der Phasengrenze als

Funktion der H¨

ohe h¨

uber der Blasenbasis. Graue Quadrate: Dynamischer Druck auf die

linke Blasenseite. Graue Rauten: Dynamischer Druck auf die rechte Blasenseite. Schwarze

Kreise: Differenz der beiden Dr¨

ucke.

13.1. Zweidimensionale Blasen

Der hydrostatische Druck verursacht die Auftriebskraft der Blase, die ausschließlich

in vertikaler Richtung wirkt. Der dynamische Druck hingegen kann auch horizontale

Kr¨

afte verursachen, wenn in dieser Richtung Asymmetrien auftreten. In den Abbil-

dungen 13.2 und 13.3 sind die Druckfelder von pdum die kleine und die große Blase

dargestellt. Wegen der Oberfl¨

achenspannung ist der Druck in der Blase hoch und

ubersteigt die Grauskala. Daher sind die Blasen in diesen Abbildungen schwarz. Un-

terdruckgebiete hingegen zeichnen sich hell ab. Sie befinden sich an der Phasengrenze

und sind asymetrisch angeordnet. Bei der kleinen Blase tritt ein Unterdruckgebiet auf

der rechten Seite auf, wo sie von der Fl¨

ussigphase st¨

arker angestr¨

omt wird. Bei der

großen Blase treten beiderseits Unterdruckgebiete auf, wovon jedoch das linke st¨

arker

ausgepr¨

agt ist. Um die Kraft des dynamischen Druckes auf die Blase zu berechnen,

muß der Beitrag jedes Fl¨

achenelementes bestimmt werden. F¨

ur ein Fl¨

achenelement

der L¨

ange sist dies F=pdsn. Die horizontale Komponente hiervon ist

Fy=pdsny=pdscos(α) = pds∆h(13.5)

mit ∆hals H¨

ohendifferenz der Endpunkte des Linienst¨

ucks s(vgl. Abbildung 13.4).

Die Horizontalkraft auf eine Blasenseite ist daher das Integral Rpd(h)dh ¨

uber die

gesamte Blasenh¨

ohe Hvon der Blasenbasis bis zu ihrem h¨

ochsten Punkt. Seien pd,L

und pd,R jeweils der Druck auf der linken und der rechten Blasenseite, so ist die

∆hs

αα

pd,L(h )

pd,R(h )

Abbildung 13.4.: Links: Oberfl¨

achenelement an der Phasengrenze einer 2D-Blase

Rechts: Phasengrenze einer 2D-Blase und Kontrollinie (gestrichelt) außerhalb der Phasen-

grenze.

13. Untersuchung der Geschwindigkeits- und Druckfelder an der Phasengrenze

horizontale Komponente des dynamischen Druckes auf die gesamte Phasengrenze

Fy=ZH

[pd,L(h)−pd,R(h)]dh . (13.6)

Weil die Oberfl¨

achenspannung in FS3D als Volumenkraft modelliert ist, sind die

Druckwerte in der Fl¨

ussigphase in der N¨

ahe der Phasengrenze stark erh¨

oht - auch in

Zellen, die keine Gasphase enthalten (f=0). Erst in einem Abstand von drei Zellbrei-

ten jenseits der Phasengrenze verschwindet der Einfluß der ”verschmierten“ Ober-

fl¨

achenspannungskraft. Daher werden die Druckwerte auf einer Kontrollinie im Ab-

stand von drei Zellen außerhalb der Phasengrenze untersucht. In Abbildung 13.2 und

13.3 sind jeweils rechts die so gewonnenen Werte f¨

ur pd,L und pd,L zusammen mit ihrer

Differenz gegen¨

uber der H¨

ohe habgetragen. In Abbildung 13.2 ist pd,L −pd,R ¨

uberwie-

gend positiv, in Abbildung 13.3 ¨

uberwiegend negativ. Es ergibt sich nach Gleichung

(13.6) in beiden F¨

allen eine Kraft in Richtung der Blasenmigration. Dies legt die Ver-

mutung nahe, daß der dynamische Druck zur Liftkraft beitr¨

agt. Wegen des starken

Wandeinflusses in 2D-Simulationen wird eine genauere quantitative Untersuchung zu

den Kr¨

aften an der Phasengrenze nur f¨

ur 3D-Blasen durchgef¨

uhrt.

13.2. Dreidimensionale Blasen

13.2.1. Zirkulation um dreidimensionale Blasen

Analog zu Abschnitt 13.1.1 soll auch f¨

ur 3D-Blasen ¨

uberpr¨

uft werden, ob die Mi-

grationsrichtung und das Vorzeichen der Zirkulation um die Blase zusammenh¨

angen.

Dazu muß die Definition der Zirkulation (13.1) von zwei auf drei Dimensionen er-

weitert werden. Nach dem Satz von Stokes l¨

aßt sich ein Linienintegral entlang der

geschlossenen Kurve Γ in ein Fl¨

achenintegral ¨

uber die von Γ umschlossene Fl¨

ache

umwandeln:

Z=IΓ

vds=ZA

(∇×v)da.(13.7)

Durch

Z=ZV

(∇×v)dv(13.8)

wird eine vektorielle Zirkulation um ein geschlossenes, endliches Volumen V(beispiels-

weise dem Blasenvolumen) definiert. Um die physikalische Bedeutung des Vektors Z

13.2. Dreidimensionale Blasen

zu verdeutlichen, soll nun eine Komponente gesondert betrachtet werden:

Zy=eyZ=ZV

(∇×v)eydv=Z ZAy

(∇×v)dady=Z IΓy

vdsdy . (13.9)

Hierbei bezeichnet Ayeine Querschnittsfl¨

ache des Volumens Vparallel zur x-z-Ebene

durch den Punkt (0,y,0). Der Rand von Ayist die geschlossenen Kurve Γy. Damit kann

Zyfolgendermaßen interpretiert werden: Sei Z(y) die Zirkulation um die Schnittfl¨

ache

Aymit dem Rand Γy, so gilt:

Zy=ZZ(y)dy . (13.10)

Die Integrationsgrenzen sind die minimale und die maximale y-Koordinate des Volu-

mens V. Wenn Vbeispielsweise eine ellipsoide Blase ist, so wird y¨

uber dem Intervall

[−dH,+dH] integriert. Weil die x-z-Ebene in der Simulation eine Symmetrieebene ist,

kann nur die y-Komponente von Zvon Null verschieden sein. Denn Schnittfl¨

achen

durch das Blasenvolumen senkrecht zur x-z-Ebene k¨

onnen wegen der Spiegelsym-

metrie des Geschwindigkeitsfeldes keine von Null verschiedene Zirkulation an ihrem

Rand aufweisen. Analog zum 2D-Fall kann eine durchschnittliche Rotation in der

Gasphase definiert werden, wenn man die Zirkulation durch das Blasenvolumen VB

teilt.

ωg=Zy

(13.11)

F¨

ur starre K¨

orper (z.B. Tragfl¨

achen), die von einem Fluid mit einer Geschwindig-

keit vrangestr¨

omt werden, ergibt sich nach dem Theorem von Kutta-Joukowski eine

Auftriebskraft FL, die proportional zur Zirkulation Zist:

FL=ρlvr×Z.(13.12)

Durch Gleichsetzen dieses Theorems mit Gleichung (11.2) ergibt sich:

CL=ωg

ωl

.(13.13)

Aus Abbildung 13.5 ist erkennbar, daß das Vorzeichen des Liftkoeffizienten nicht von

dem Vorzeichen von ωgabh¨

angt. Insbesondere ist Gleichung (13.13) nicht erf¨

ullt. Es

l¨

aßt sich auch kein stoffunabh¨

angiger funktionaler Zusammenhang zwischen CLund

ωgherstellen.

13. Untersuchung der Geschwindigkeits- und Druckfelder an der Phasengrenze

-1,4 - 1,2 -1,0 - 0,8 -0,6 -0,4 -0,2 0,0 0,2 0,4 0,6 0,8

-0,6

-0,4

-0,2

0,0

0,2

0,4

19 mPa s 53 mPa s

y = x

Abbildung 13.5.: Liftkoeffizient in Abh¨

angigkeit von ωg/ωlf¨

ur Luftblasen in Glyce-

rinl¨

osungen (sie auch Tabelle 12.1, sowie Abbildung 12.7).

1.2 1.4 1.6 1.8

0.8

1.0

1.2

1.4

1.6

−10.0

−6.7

−3.3

0.0

3.3

6.7

10.0

13.3

16.7

20.0

dyn. pres.

[Pa]

t = 0.200000 s

cycle = 1

1.5 2.0 2.5

1.5

2.0

2.5

3.0

−21.7

−14.4

−7.2

0.0

7.2

14.4

21.7

28.9

36.1

dyn. pres.

[Pa]

t = 0.200000 s

cycle = 1

Abbildung 13.6.: Dynamischer Druck um 3D-Blasen, µl= 53 mPa s. Links: Kleine Blase,

d=3,2 mm. Rechts: Große Blase, d=5,6 mm. Scherrate ω= 8.3s−1. Die Gegenstr¨

omung

nimmt jeweils nach rechts zu. Dazugeh¨

orige Geschwindigkeitsfelder: Siehe Abbildung 12.9.

100

13.2. Dreidimensionale Blasen

13.2.2. Dynamischer Druck und Scherkraft an der Phasengrenze

Dreidimensionale Blasen sind entlang ihres ”¨

Aquators“ von einem ringf¨

ormigen Un-

terdruckgebiet umgeben, welches sie in horizontaler Richtung dehnt (Abbildung 13.6).

Das verdr¨

angte Wasser streicht an der Blase vorbei, wobei die Wassergeschwindigkeit

zunimmt. Gem¨

aß dem Bernoulli-Effekt sinkt der Wasserdruck mit dem Geschwindig-

keitsanstieg. Dies erkl¨

art auch, warum schnelle Blasen flacher sind, als langsame (vgl.

Abb. 12.2 und 12.3). F¨

ur alle 3D-Blasen bei µl≤53 mPa s liegt das absolute Mini-

mum des dynamischen Druckes auf der st¨

arker angestr¨

omten Seite, unabh¨

angig vom

Blasendurchmesser. Nur f¨

ur große Blasen bei h¨

oheren Viskosit¨

aten (z.B. d=10 mm

und µl= 100 mPa s) tritt das absolute Minimum auf der schw¨

acher angestr¨

omten

Seite auf. Um die gesamte Kraft Fpauszurechnen, die der dynamische Druck pdauf

die Blase aus¨

ubt, muß pd¨

uber die Phasengrenze Γ integriert werden:

Fp=ZΓ−pdndA , (13.14)

mit der ¨

außeren Fl¨

achennormalen n. Wie bereits in Abschnitt 13.1.2 dargelegt, sind

die durch FS3D errechneten Druckwerte in der N¨

ahe der Phasengrenze durch die

Modellierung der Oberfl¨

achenspannung als Volumenkraft stark erh¨

oht. Erst in einem

Abstand von drei bis vier Zellbreiten von der Phasengrenze verschwindet der Einfluß

dieser Volumenkraft, und die errechneten Werte k¨

onnen als physikalische Druckwerte

in der Fl¨

ussigphase interpretiert werden. Daher kann Fpnicht direkt durch Integrati-

on ¨

uber die Phasengrenze berechnet werden, sondern wird durch eine Extrapolation

angen¨

ahert. Hierzu werden Fl¨

achen Aδermittelt, die in einem festen Abstand δpar-

allel zur Phasengrenze in der Fl¨

ussigphase verlaufen (vgl. Abbildung 13.7). Der Algo-

rithmus zur Ermittlung dieser Fl¨

achen ist im Anhang A.3 beschrieben. Im Grenzfall

Aδ

Abbildung 13.7.: Fl¨

ache Aδim Abstand δzur Phasengrenze Γ.

101

13. Untersuchung der Geschwindigkeits- und Druckfelder an der Phasengrenze

δ→0 geht Aδin die Phasengrenze Γ ¨

uber. Daher gilt

Fp= lim

δ→0Fp,δ mit Fp,δ =ZAδ−pdndA . (13.15)

Analog zum Druck kann auch die Scherkraft FSan der Phasengrenze untersucht

werden:

FS=ZΓ

S n dAmit S=µ(∇v+∇vT).(13.16)

Anders als beim Druck weist das Geschwindigkeitsfeld keine Unstetigkeit an der Pha-

sengrenze auf. Allerdings ergibt sich durch den Sprung der Viskosit¨

at an der Phasen-

grenze eine Unstetigkeit f¨

ur ∇v. Daher wird auch f¨

ur FSeine Extrapolation durch-

gef¨

uhrt:

FS= lim

δ→0FS,δ mit FS,δ =ZAδ

S n dA . (13.17)

13.2.3. Validierung des Verfahrens zur Berechnung der Kr¨

afte an der

Phasengrenze

Um das numerische Verfahren zur Berechnung von Fpund FSzu validieren, werden

zun¨

achst Blasen in ruhenden Fl¨

ussigkeiten (d.h. bei ω= 0) untersucht, die nach einer

anf¨

anglichen Beschleunigungsphase mit konstanter Geschwindigkeit wTaufsteigen.

Nach Erreichen der Endgeschwindigkeit wTsind Widerstandskraft FDund Auftriebs-

kraft FAim Gleichgewicht:

FD=−FA= ∆ρVBg.

Die Widerstandskraft setzt sich aus den Beitr¨

agen des dynamischen Druckes und der

Scherkraft zusammen: FD=Fp+FS. Es muß also f¨

ur jede Blase nach Erreichen der

Endgeschwindigkeit wTgelten:

Fp+FS=−FA.(13.18)

Diese Gleichheit wird f¨

ur verschiedene Blasendurchmesser dbei µl= 19 mPa s

und µl= 53 mPa s ¨

uberpr¨

uft. In Abbildung 13.8 wird am Beispiel einer Blase mit

d=3,2 mm bei µl= 53 dargestellt, wie Fpund FSdurch Extrapolation gewonnen wer-

den. Als Abstand zur Phasengrenze werden jeweils Vielfache der Zellbreite cgew¨

ahlt:

102

13.2. Dreidimensionale Blasen

0 2 4 6 8 10

0,2

0,3

0,4

0,5

0,6

0,7

|F | / |F

Abbildung 13.8.: Druckkraft und Scherkraft auf Kontrollfl¨

achen Aδin Abh¨

angigkeit vom

Abstand δzur Phasengrenze. Die Kr¨

afte sind durch die Auftriebskraft normiert, der Ab-

stand ist durch die Kantenl¨

ange cder kubischen Zellen normiert.

δ= 9c, 8c, ..., 0. Die Werte von Fp,δ und FS,δ lassen sich f¨

ur δ > 3cjeweils sehr

gut durch ein quadratisches Polynom fitten. Extrapoliert man die Fits in Richtung

δ→0 so ergibt sich in diesem Falle Fp=−0,36 FAund FS=−0,62 FAund damit

insgesamt FD=−0,98 FA. F¨

ur alle simulierten Blasen ist Gleichung (13.18) (bis

auf relative Fehler von wenigen Prozent) erf¨

ullt, wenn beide Anteile des Drags durch

Extrapolation bestimmt werden. Berechnet man Fpoder FSoder beide Kr¨

afte statt-

dessen direkt an der Phasengrenze (δ= 0), so ist Gleichung (13.18) im Allgemeinen

nicht erf¨

ullt.

Aus Abbildung 13.9 ist ersichtlich, daß im Falle Re→0 der Anteil von Fpan der ge-

samten Widerstandsraft FDein Drittel betr¨

agt. Dieser Wert ergibt sich auch aus der

analytischen L¨

osung von Hadamard (1911) und Rybczynski (1911) f¨

ur kugelf¨

ormige

Fluidpartikel im Grenzfall µg→0. Je gr¨

oßer hingegen die Reynoldszahl wird, desto

gr¨

oßer ist der Anteil des dynamischen Druckes an der Widerstandskraft. Auch dies

entspricht der Theorie, denn im Grenzfall Re→ ∞ spielen Scherkr¨

afte keine Rolle.

Durch Bilanzierung der Kr¨

afte an frei aufsteigenden Blasen konnten die numerischen

103

13. Untersuchung der Geschwindigkeits- und Druckfelder an der Phasengrenze

Routinen zur Berechnung der Oberfl¨

achenintegrale validiert werden. Es wurde gezeigt,

daß die durch Extrapolation ermittelten Werte f¨

ur die Kr¨

afte an der Phasengrenze

plausibel sind.

0 10 20 30 40 50 60 70 80

0,3

0,4

0,5

0,6

0,7

0,8

53 mPa s

19 mPa s

| / | F

Abbildung 13.9.: Anteil des dynamischen Druckes an der Widerstandskraft in Abh¨

angig-

keit von der Reynoldszahl f¨

ur Blasen verschiedener Gr¨

oßen bei µl= 19 und 53 mPa s.

Berechnung der Horizontalkr¨

afte an Blasenoberfl¨

achen in Scherstr¨

omungen

Betrachtet man die Druckfelder in Abbildung 13.6, so liegt die Vermutung nahe, daß

sich im 3D-Fall aus dem dynamischen Druck allein nicht die unterschiedliche Migra-

tionsrichtung kleiner und großer Blasen in Scherstr¨

omungen erkl¨

aren l¨

aßt. Allerdings

kann man die horizontale Komponente der Druckkraft Fpschlecht per Augenschein

absch¨

atzen, weil die Seiten der großen Blase unterschiedlich gekr¨

ummt sind. Daher

werden die horizontalen Komponenten von Fpund FSmit dem oben beschriebenen

Verfahren berechnet. Hierf¨

ur werden die Geschwindigkeits- und Druckfelder ausge-

wertet, die sich jeweils nach dem Erreichen konstanter Blasen-Relativgeschwindigkeit

einstellen. Weil die horizontale Blasengeschwindigkeit dann konstant ist, m¨

ussen sich

104

13.2. Dreidimensionale Blasen

die Horizontalkr¨

afte an der Phasengrenze im Gleichgewicht befinden:

FS,x =−Fp,x .(13.19)

In den Abbildungen 13.10 und 13.11 sind die Horizontalkr¨

afte in Abh¨

angigkeit vom

Blasendurchmesser aufgetragen. Bei µl=53 mPa s ist Gleichung (13.19) bis auf Feh-

ler von wenigen Prozent erf¨

ullt. Fp,x ist stets positiv, d.h. große Blasen migrieren

entgegengesetzt zu der von Fp,x vorgegebenen Richtung. Bei µl=19 mPa s stimmt

die Migrationsrichtung stets mit der Richtung von Fp,x ¨

uberein, jedoch ist Gleichung

(13.19) nicht erf¨

ullt. Zwar ist das Vorzeichen von FS,x stets dem von Fp,x entgegen-

gesetzt, jedoch ist der Betrag von FS,x nur etwa halb so groß wie der von Fp,x. Die

Blase m¨

ußte also in horizontaler Richtung beschleunigen, was aber nicht der Fall ist.

M¨

oglicherweise versagt die Extrapolation hier, weil die Blasen bei µl=19 mPa s sehr

flach sind, und die Integration um stark gekr¨

ummte Oberfl¨

achen aufgrund der be-

schr¨

ankten Gitteraufl¨

osung ungenau wird. Insgesamt ergibt sich aus den Werten f¨

die Horizontalkr¨

afte an der Phasengrenze kein schl¨

ussiges Gesamtbild. Lediglich f¨

Blasen mit positivem Liftkoeffizienten kann man feststellen, daß sie stets in die Rich-

tung migrieren, die aus dem dynamische Druck auf der Phasengrenze resultiert. Sie

migrieren zur st¨

arker angestr¨

omten Seite, wo wegen des Bernoulli-Effekts der Unter-

druck am st¨

arksten ist. Das Vorzeichen der Scherkraft FS,x ist f¨

ur alle Blasengr¨

oßen

dem der Druckkraft Fp,x entgegengesetzt.

105

13. Untersuchung der Geschwindigkeits- und Druckfelder an der Phasengrenze

-0,4 -0,3 -0,2 -0,1 0,0 0,1 0,2 0,3

-0,03

-0,02

-0,01

0,00

0,01

0,02

0,03

S,x

p,x

/ |F

Abbildung 13.10.: Horizontale Komponenten der Druck- und Scherkraft auf Blasen bei

µl= 53 mPa s in Abh¨

angigkeit vom Liftkoeffizienten. Die Horizontalkr¨

afte sind durch die

jeweilige Auftriebskraft normiert.

-0,7 -0,6 -0,5 -0,4 -0,3 -0,2 -0,1 0,0 0,1 0,2 0,3 0,4

-0,04

-0,03

-0,02

-0,01

0,00

0,01

0,02

0,03

0,04

0,05

p,x

S,x

/ |F

Abbildung 13.11.: Horizontale Komponenten der Druck- und Scherkraft auf Blasen bei

µl= 19 mPa s in Abh¨

angigkeit vom Liftkoeffizienten. Die Horizontalkr¨

afte sind durch die

jeweilige Auftriebskraft normiert.

106

14. Horizontalgeschwindigkeiten im Blasennachlauf

In Abbildung 12.10 und 13.1 ist zu erkennen, daß der Nachlauf großer Blasen leicht in

Richtung st¨

arkerer Gegenstr¨

omung abgeschr¨

agt ist. Tomiyama (1992) vermutet, daß

die Geometrie des Nachlaufes wesentlich zum Lift beitr¨

agt. In einem abgeschr¨

agten

Nachlauf findet ein horizontaler Fl¨

ussigkeitstransport statt. Eingehende numerische

Untersuchungen zur Struktur des Nachlaufes hinter Partikeln in linearen Scherstr¨

omun-

gen wurden von Richard et al. (2007) durchgef¨

uhrt. Als Modell f¨

ur Blasen wurden

formfeste, abgeflachte Ellipsoide mit einer Slip-Randbedingung verwendet. Die rotati-

onssymmetrischen Elipsoide wurden fest im Rechengebiet positioniert, und die Kr¨

afte

der umgebenden Scherstr¨

omung auf den K¨

orper in Abh¨

angigkeit von der Reynolds-

zahl, der Scherrate und des Durchmesserverh¨

altnisses berechnet. Aufgrund dieser Vor-

gehensweise sind die Resultate dieser Rechnungen nicht direkt mit experimentellen

Ergebnissen zu vergleichen, weil das Durchmesserverh¨

altnis und die Reynoldszahl im

Experiment nicht unabh¨

angig voneinander variiert werden k¨

onnen. Blasenform und

Geschwindigkeit sind im Experiment wie auch in Simulationen mit freier Phasengren-

ze untrennbar verkn¨

upft: Je nach Blasengr¨

oße und Scherrate stellt sich eine Form

und eine Aufstiegsgeschwindigkeit ein. Daher werden im Folgenden die Ergebnisse

von Richard et al. (2007) vorwiegend qualitativ mit den experimentellen Resultaten

von Tomiyama und der eigenen Rechnungen mit FS3D verglichen. Die Liftkraft ist

in den Rechnungen von Richard et al. (2007) abh¨

angig von der Scherrate. Negative

Liftkr¨

afte treten generell erst bei Durchmesserverh¨

altnissen von dV/dH<0.455 auf.

Je h¨

oher die Scherrate ist, desto kleiner ist das kritische Durchmesserverh¨

altnis, bei

dem die Liftkraft ihr Vorzeichen wechselt. Im Experiment von Tomiyama wie auch in

den FS3D-Simulationen liegt das kritische Durchmesserverh¨

altnis zwischen 0,55 und

0,65 - je nach Stoffsystem (vgl. Abbildung 11.5 und 12.13). Eine Abh¨

angigkeit der

Liftkraft von der Scherrate besteht hier nicht (vgl. Abschnitt 12.5). Das Vorzeichen

der Liftkraft bringen Richard et al. (2007) mit der Rotation in der Fl¨

ussigphase hinter

der Blase in Zusammenhang. F¨

ur ein Ellipsoid mit einem Durchmesserverh¨

altnis von

107

14. Horizontalgeschwindigkeiten im Blasennachlauf

1.00 1.10 1.20 1.30 1.40 1.50 1.60

1.00

1.10

1.20

1.30

1.40

vof−function

[−]

vmin = 0.551725

vmax = 23.3500

velocity

[cm/s]

t = 0.100001 s

c = 3081

2.0 2.5 3.0 3.5

1.6

1.8

2.0

2.2

2.4

2.6

2.8

Abbildung 14.1.: Geschwindigkeitsfelder auf der Symmetrieebene f¨

ur Blasen in einer Gly-

cerinl¨

osung (µl= 19 mPa s). Die Gegenstr¨

omung nimmt nach rechts zu. Links: Kleine Blase

(d=3,2 mm), Migration nach rechts. Rechts: Große Blase (d=5,6 mm), Migration nach links.

dV/dH= 0.444 beobachten sie eine von Null verschiedene Rotation ωz=ez·rotvlim

Nachlauf. Aus Symmetriegr¨

unden ist das Vorzeichen der Rotation im Halbraum y < 0

stets umgekehrt wie im Halbraum y > 0. Daraus ergibt sich ein Transport der Fl¨

ussig-

phase entlang der Symmetrieebene (hier die x-z-Ebene). Gleichzeitig ist eine Liftkraft

in entgegengesetzte Richtung zu beobachten. Abh¨

angig von der Reynoldszahl kehrt

sich das Vorzeichen der Rotation ωzum, und mit ihr zugleich auch das Vorzeichen des

Liftkoeffizienten. Diese Beobachtung l¨

aßt sich zu folgender anschaulichen Erkl¨

arung

zusammenfassen: Hinter der Blase entsteht ein Paar von entgengesetzt rotierenden

Wirbelwalzen, zwischen denen Fl¨

ussigkeit in Richtung ±extransportiert wird. Die

Blase migriert stets in entgegengesetzte Richtung, als sei sie von der transportier-

ten Fl¨

ussigkeit wie von einem Jet angetrieben. Auch in analytischen Rechnungen

zu sph¨

arischen Partikeln in Scherstr¨

omungen reibungsfreier Fluide ergibt sich eine

Rotation ωzim Nachlauf (Auton, 1984; Legrende und Magnaudet, 1998). Unter Ver-

wendung der Wirbeltransportgleichung bei Vernachl¨

assigung des Terms f¨

ur viskose

Diffusion leiten Richard et al. (2007) die Entstehung der Rotation im Nachlauf sph¨

ari-

scher Partikel her und motivieren, wie eine elliptische Deformation des Partikels zur

Umkehrung der Rotation f¨

uhren kann. Diese ¨

Uberlegungen sind wegen der Annahme

108

Draufsicht

−

Seitansicht

a) b)

Abbildung 14.2.: Rotation ωzim Nachlauf einer Blase. Gestrichelte Linie: Symmetrieebe-

ne (x-z Ebene). a) Vertikaler Schnitt durch den Blasenschwerpunkt (y-z Ebene).

b) Horizontaler Schnitt durch den Nachlauf. ¨

Außerer Kreis: Projektion des Blasenrandes auf

die horizontale Ebene. Innere Kreise: Gegensinnig rotierende Wirbelwalzen (schematisch).

der Reibungsfreiheit nur im Fall Re≫1 g¨

ultig. Auch die numerischen Untersuchungen

von Richard et al. (2007) konzentrieren sich ¨

uberwiegend auf den Bereich Re≥400.

F¨

ur die in der vorliegenden Arbeit untersuchten Blasen in Glycerinl¨

osungen sind die

Reynoldszahlen jedoch relativ klein (Re<100) und die Formen weichen teilweise

erkennbar von Ellipsoiden ab. Im Folgenden soll ¨

uberpr¨

uft werden, ob hier auch die

oben beschriebenen Wirbelwalzen im Nachlauf auftreten und ob ihr Drehsinn mit dem

Vorzeichen des Liftkoeffizienten in Verbindung steht. In Abbildung 14.3 sind f¨

ur eine

Blase mit CL<0 (d= 5,6 mm, µl= 19 mPa s) die Rotation ωzund die Horizontal-

geschwindigeit uauf der x-z-Ebene dargestellt. Man beachte, daß diese Schnittebene

senkrecht zur Symmetrieebene ist. Die Druck- und Geschwindigkeitsfelder wurden

in den vorangegangenen Kapiteln stets ausschließlich auf der Symmetrieebene (x-z-

Ebene) untersucht. Das Vorzeichen von ωzist im Nachlauf positiv, wodurch auch u

im Nachlauf positiv ist. Die Blase migriert in entgegengesetzte Richtung, in ¨

Uber-

einstimmung mit dem von Richard et al. (2007) beschriebenen Mechanismus. Anders

als in der schematischen Zeichnung 14.2b liegt das Zentrum der Wirbelwalze in x-

Richtung nicht mittig unter der Blase, sondern in Richtung st¨

arkerer Gegenstr¨

omung

verschoben.

109

14. Horizontalgeschwindigkeiten im Blasennachlauf

Aus Abbildung 14.5 ist ersichtlich, daß die Horizontalgeschwindigkeit uim Nachlauf

der großen Blase (d= 5,6 mm) fast ¨

uberall positiv ist - mit Ausnahme der unmit-

telbaren N¨

ahe der Phasengrenze. Damit wird im Nachlauf Fl¨

ussigkeit vorwiegend in

Richtung extransportiert. Die große Blase (d= 5,6 mm) migriert in entgegengesetzte

Richtung in ¨

Ubereinstimmung mit der Beobachtung von Richard et al. (2007). Anders

verh¨

alt es sich bei der kleinen Blase (d= 2,8 mm), die in der gleichen Fl¨

ussigkeit

(µl= 19 mPa s) einen positiven Liftkoeffizienten hat. In ihrem Nachlauf ist u¨

uberall

positiv (vgl. Abbildung 14.5). Demnach wird im Nachlauf Fl¨

ussigkeit in die gleiche

Richtung transportiert, in die auch die Blase migriert. Das Erkl¨

arungsmodell f¨

ur die

Vorzeichenumkehr der Liftkraft von Richard et al. (2007) trifft auf dieses Blasenpaar

also nicht zu.

1.5 2.0 2.5 3.0 3.5 4.0

0.0

0.2

0.4

0.6

0.8

−8.2

−5.5

−2.7

0.0

2.7

5.5

8.2

10.9

13.6

16.4

19.1

[1/s]

t = 0.300000 s

cycle = 1

+20

+10

−10

1.5 2.0 2.5 3.0 3.5 4.0

0.0

0.2

0.4

0.6

0.8

−2.9

−2.5

−2.1

−1.7

−1.2

−0.8

−0.4

0.0

0.4

0.8

1.2

[cm/s]

t = 0.300000 s

−3

−1

−2

Abbildung 14.3.: Vertikaler Schnitt durch den Schwerpunkt einer Blase (d=5,6 mm, µl=

19 mPa s). Elliptische Line: Blasenrand. Es ist nur die Halbebene y > 0 dargestellt. Links:

Rotation ωz(vgl. Abb. 14.2 a) Rechts: Horizontale Geschwindigkeiskomponente vx.

110

1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0

0.0

0.2

0.4

0.6

0.8

1.0

1.2

0.0

1.1

2.2

3.3

4.4

5.5

6.6

7.8

8.9

10.0

11.1

[1/s]

t = 0.300000 s

cycle = 1

1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0

0.0

0.2

0.4

0.6

0.8

1.0

Abbildung 14.4.: Horizontaler Schnitt durch den Nachlauf einer Blase (d=5,6mm, µl=

19 mPa s). Die Schnittebene liegt 1,5dunterhalb des Blasenschwerpunkts. Es ist nur

die Halbebene y > 0 dargestellt. Kreisf¨

ormige Line: Projektion des Blasenrandes auf die

Schnittebene. Links: Rotation ωz(vgl. Abb. 14.2 b) Rechts: Horizontale Geschwindigkeiten.

-15 -10 -5 0 5 10 15

-2,5

-2,0

-1,5

-1,0

-0,5

0,0

0,5

1,0

1,5

[cm/ s]

-2 z / d

d=2,8 mm

d=5,6mm

Nachlauf

Vorlauf

Abbildung 14.5.: Horizontalgeschwindigkeit uin Abh¨

angigkeit von der H¨

ohe zauf einer

vertikalen Linie durch den Blasenschwerpunkt. Die H¨

ohe wird unter Verwendung von dVso

normiert, daß die Phasengrenze f¨

ur beide Blasen (d=3,2 mm und d=5,6 mm) bei ±1 liegt.

111

14. Horizontalgeschwindigkeiten im Blasennachlauf

112

15. Zusammenfassung des zweiten Teils

Neben experimentellen Untersuchungen, werden h¨

aufig auch numerische Str¨

omungs-

simulationen eingesetzt, um das Verst¨

andnis von mehrphasigen Str¨

omungen in tech-

nischen Anlagen - wie beispielsweise in Blasens¨

aulen - zu verbessern. Eine akkura-

te Simulation der Blasenbewegung erfordert eine detaillierte Kenntnis der Wechsel-

wirkungskr¨

afte zwischen Gas- und Fl¨

ussigphase. Steigen Blasen in nicht homogenen

Str¨

omungen auf, so tritt eine Kraft senkrecht zu ihrer Relativgeschwindigkeit ge-

gen¨

uber der Fl¨

ussigphase auf. Diese Kraft wird Liftkraft genannt. Mit dem Zweipha-

sencode FS3D wird die Bewegung von Blasen in linearen Scherstr¨

omungen simuliert.

Aus der Analyse der Bewegungsbahnen wird die jeweils aufgetretene Liftkraft be-

stimmt, deren Vorzeichen und Betrag in Form eines dimensionslosen Liftkoeffizienten

beschrieben wird. Die Abh¨

angigkeit des Liftkoeffizienten von der Blasengr¨

oße und den

Stoffeigenschaften der Fl¨

ussigkeit wird numerisch untersucht und mit experimentellen

Resultaten verglichen (Tomiyama et al., 2002). Wie im Experiment h¨

angt auch in den

Simulationen der Liftkoeffizent im Wesentlichen von einer modifizierten E¨

otv¨

os-Zahl

ab, die auf dem horizontalen Durchmesser der Blase basiert. Die E¨

otv¨

os-Zahl ist eine

dimensionslose Kennzahl, die das Verh¨

altnis von Auftriebskraft zu Oberfl¨

achenspan-

nung beschreibt. Blasen mit einer kleinen E¨

otv¨

os-Zahl sind eher kugelf¨

ormig, wohinge-

gen Blasen mit einer großen E¨

otv¨

os-Zahl st¨

arker abgeflacht sind. In ¨

Ubereinstimmung

mit dem Experiment ist der Liftkoeffizient f¨

ur Blasen mit kleiner E¨

otv¨

os-Zahl positiv

und f¨

ur Blasen mit großer E¨

otv¨

os-Zahl negativ. Dies bedeutet eine Richtungsumkehr

der Liftkraft. In Experimenten mit Blasen in Glycerinl¨

osungen variierten Tomiyama

et al. (2002) die Blasengr¨

oßen. Diese Experimente wurden nachgerechnet. Zus¨

atzlich

wurde in den Simulationen Oberfl¨

achenspannung und Gravitation variiert. F¨

ur alle

F¨

alle konnte gezeigt werden, daß sich der Liftkoeffizient CLals eine Funktion der

modifizierten E¨

otv¨

oszahl E¨

oHbeschreiben l¨

aßt. Die ermittelten Funktionen CL(E¨

oH)

sind sich ¨

ahnlich, unabh¨

angig davon, auf welche Weise E¨

oHvariiert wurde. Die aus

den Simulationen ermittelten Liftkoeffizienten sind jedoch stets etwas kleiner als die

113

15. Zusammenfassung des zweiten Teils

experimentell bestimmten. Diese systematische Abweichung wurde auch in anderen

numerischen Rechnungen beobachtet (Rusche, 2002; Kuipers et al., 2006; Sousa et al.,

2006). M¨

oglichweise treten im Experiment Ph¨

anomene auf, die in den numerischen

Codes nicht modelliert sind. Beispielsweise spielen im Experiment Verunreinigungen

der Fl¨

ussigkeit eine Rolle, welche sich an der Blasenoberfl¨

ache anlagern und die Ober-

fl¨

acheneigenschaften ver¨

andern (vgl. Kapitel 10).

Die Simulationen liefern vollst¨

andige Geschwindigkeits- und Druckfelder. Durch Aus-

wertung dieser Felder kann gezeigt werden, daß die Liftkraft an dreidimensionalen

Blasen nicht als Funktion der Zirkulation des Geschwindigkeitsfeldes an der Phasen-

grenze beschrieben werden kann. Demnach kann die Liftkraft nicht auf den Magnus-

effekt zur¨

uckgef¨

uhrt werden, der an umstr¨

omten rotierenden Zylindern auftritt.

Es wurden Routinen zur Integration und Extrapolation der Druck- und Scherkr¨

afte

an der Phasengrenze entwickelt und f¨

ur frei aufsteigende Blasen validiert. F¨

ur eher

kugelf¨

ormige Blasen (E¨

oHklein) in linearen Scherstr¨

omungen kann auf diese Weise

gezeigt werden, daß der dynamische Druck um die Blasen in Richtung der Liftkraft

wirkt. Dies l¨

aßt sich durch den Bernoulli-Effekt erkl¨

aren. F¨

ur st¨

arker deformierte Bla-

sen (E¨

oHgroß) wird kein Zusammenhang zwischen der Richtung der Liftkraft und

dem dynamischen Druck in der Umgebung der Blase festgestellt. Insgesamt kann mit

den verwendeten Methoden nicht eindeutig gezeigt werden, daß die Liftkraft in allen

F¨

allen durch den dynamischen Druck um die Blase verursacht wird. Weil die Phasen-

grenze in der VoF-Methode nicht r¨

aumlich scharf lokalisiert ist, sondern ¨

uber einige

Zellbreiten verschmiert wird, sind exakte Werte f¨

ur Druck und Geschwindigkeit an

der Phasengrenze nicht eindeutig zu ermitteln. Daher bietet es sich an, f¨

ur Untersu-

chungen der Kr¨

afte an der Phasengrenze einen numerischen L¨

oser zu verwenden, in

dem die Phasengrenze scharf lokalisiert ist - beispielsweise durch ein bewegtes Gitter.

114

A. Anhang

A.1. Symbolverzeichnis

Tabelle A.1.: Lateinische Symbole

Symbol Einheit Bedeutung

CD[ - ] Widerstandsbeiwert (drag coefficient)

dm¨

Aquivalenzdurchmesser

dWm Wandabstand

f[-] Volumenanteilsfunktion (”VoF-Variable“)

gm/s2Erdbeschleunigung

gm/s2Erdbeschleunigung (vektoriell)

I[-] Einheitstensor

nn Fl¨

achennormale

pPa Druck

pdPa dynamischer Druck

S[-] Z¨

ahigkeitsspannungstensor

T[-] Spannungstensor

um/s horizontale Geschwindigkeit (x-Komponente)

um/s Geschwindigkeitsvektor

vm/s horizontale Geschwindigkeit (y-Komponente)

vm/s Geschwindigkeitsvektor

Vm3Volumen

VBm3Blasenvolumen

wm/s vertikale Geschwindigkeit (z-Komponente)

wTm/s vertikale Endgeschwindigkeit (terminal velocity)

Zm2/s Zirkulation

115

A. Anhang

Tabelle A.2.: Dimensionslose Kennzahlen

E¨

o E¨

otv¨

oszahl

Eu Eulerzahl

Fr Froudezahl

Mo Mortonzahl

Re Reynoldszahl

Sr dimensionslose Scherrate

λViskosit¨

atenverh¨

altnis

Tabelle A.3.: Griechische Symbole

Symbol Einheit Bedeutung

κ1/m2Kr¨

ummung

λKm Kolmogorov-L¨

ange

µPa s dynamische Viskosit¨

νm2/s kinematische Viskosit¨

ω1/s Scherrate

ρkg/m3Dichte

σN/m Grenzfl¨

achenspannung

Tabelle A.4.: Indices

Symbol Bedeutung

c kontinuierlich

d dispers

g gasf¨

ormig

l fl¨

ussig (liquid)

B Blase

P Partikel

H horizontal

V vertikal

r relativ

Γ Phasengrenze

116

A.2. Das Druckfeld in einem Hill’schen Wirbel

Das Geschwindigkeitsfeld eines Hill’schen Wirbels lautet in Kugelkoordinaten:

ur= (a+br2) cos θ , (A.1)

uθ=−(a+ 2br2) sin θ . (A.2)

Durch Vergleich mit (6.4) und (6.5) ergibt sich, daß das Geschwindigkeitsfeld der

Hadamard-Rybzinski-L¨

osung innerhalb des fluiden Partikels ein Hill’scher Wirbel mit

a=wT

2 + 2λund b=−a

R2(A.3)

ist. Das dazugeh¨

orige Druckfeld ergibt sich, wenn man die Impulsgleichung nach ∇p

aufl¨

ost.

∇p=−ρ∂tu−ρu·∇u+µ∆u−ρgez(A.4)

Da der Hill’sche Wirbel zeitunabh¨

angig ist, entf¨

allt der erste Term auf der rechten

Seite. Wird zus¨

atzlich eine schleichende Str¨

omung angenommen, so entf¨

allt auch der

n¨

achste Term, der den konvektiven Impulstransport beschreibt. Es verbleibt

∇p=µ∆u−ρgez.(A.5)

Nun soll ∆uberechnet werden. Dies bereitet in Kugelkoordinaten Schwierigkeiten,

weil dann der Ursprung ein ausgezeichneter Punkt ist: In Kugelkoordinaten ist ei-

ne negative radiale Koordinate rnicht zul¨

assig und somit liegt der Ursprung mit

r= 0 am Rand des Definitionsgebietes. Dieses Problem wird umgangen, indem das

Geschwindigkeitsfeld des Hill’schen Wirbels in kartesische Koordinaten umgewandelt

wird:

u=





−bxz

−byz

2bx2+ 2by2+bz2+a







.(A.6)

Die Anwendung des kartesischen Laplace-Operators ∆ := ∂2

x+∂2

y+∂2

zliefert

∆u=





10 b







.(A.7)

117

A. Anhang

Einsetzen dieses Ergebnisses in (A.5) liefert

∂xp= 0,(A.8)

∂yp= 0,(A.9)

∂zp= 10µb −ρg . (A.10)

Dies bedeutet, daß der Druck im Hill’schen Wirbel in horizontaler Richtung konstant

ist. In vertikaler Richtung w¨

achst der Druck linear mit der H¨

ohe z. Integration nach

zliefert

p(z) = 10µb z −ρg z +c=−5µ wT

λ+ 1

R2−ρgz +c . (A.11)

Wenn sich ein Hill’scher Wirbel in einem fluiden Partikel mit einer Oberfl¨

achen-

spannung σbefindet, und ein Außendruck p0auf der H¨

ohe z= 0 herrscht, so gilt

c=p0+ 2σ/R2. Durch Umwandlung in Kugelkoordinaten (z=rcos θ) geht Glei-

chung (A.11) schließlich in (6.8) ¨

uber.

Setzt man keine schleichende Str¨

omung voraus, so muß der konvektive u·∇uTerm

ber¨

ucksichtigt werden. Es ist nun ein Druck pzu finden, der folgende Gleichung erf¨

ullt:

∇p=ρu·∇u+µ∆u−ρgez.(A.12)

Diese Gleichung ist gel¨

ost, wenn man ˜pfindet mit

∇˜p=ρu·∇u.(A.13)

Dieser Druck ist dann zur bisherigen L¨

osung psnach Gleichung (A.11) zu addie-

ren. Die Summe p=ps+ ˜pist L¨

osung der Gleichung (A.12). Durch Einsetzen des

Geschwindigkeitsfeldes (A.6) in den konvektiven Term erh¨

alt man

u·∇u=





−2b2x3−2b2xy2−bxa

−2b2y3−2b2x2y−bya

2z3b2+ 2abz







.(A.14)

Durch Einsetzen dieses Ausdrucks in (A.13) und Integration ergibt sich

˜p=ρ

2b2(z4−x4−y4) + ρ

2ab(2z2−x2−y2)−ρ b2x2y2.(A.15)

118

A.2. Das Druckfeld in einem Hill’schen Wirbel

Dieses Druckfeld wird nun auf der vertikalen Linie (0,0,z) untersucht. Mit b=−a/R2

(vgl. (A.3)) ergibt sich

˜p(z) = −ρa2

2R4z4−ρa2

R2z2.(A.16)

Wenn der Hill’sche Wirbel sich in einem Tropfen mit dem Radius Rbefindet, so

erreicht der Betrag von ˜p(z) sein Maximum bei z=±R. An diesen Stellen gilt

˜p=ρ a2/2. Die Gr¨

oße dieses Beitrages des konvektiven Terms soll nun mit dem

hydrodynamischen Anteil pdaus (A.11) verglichen werden. Dieser hat an den Stellen

z=±Rden Betrag pd= 10µbR = 10µa/R. Nun wird ˜phierzu ins Verh¨

altnis gesetzt:

˜p

=ρRa

20µ=ρRwT

40(1 + λ).(A.17)

Interpretiert man ρRwT/µ als Reynoldszahl des Hill’schen Wirbels, so gilt

Re ≪40 ⇒˜p≪pd.(A.18)

Diese Bedingung ist f¨

ur alle zur Validierung des Codes simulierten Tropfen erf¨

ullt.

Deshalb tauchen im inneren Druckfeld praktisch keine Abweichungen von dem linea-

ren Druckverlauf gem¨

aß (A.11) auf.

119

A. Anhang

A.3. Konstruktion einer einh¨

ullenden Fl¨

ache parallel zur

Phasengrenze

Es liege eine Verteilung der dispersen Phase mit der Phasengrenze Γ vor. Gesucht ist

eine einh¨

ullende Fl¨

ache Γ′im Abstand δzu Γ. Zur Veranschaulichung wird der Al-

gorithmus zur L¨

osung dieses Problems an einem eindimensionalen Fallbeispiel vorge-

stellt: Die Verteilung der VoF Variable sei als eindimensionale Datenreihe figegeben,

wobei idie Zellnummer ist. In Abbildung A.1 ist ein Beispiel f¨

ur fidargestellt mit

fi>0 in den Intervallen 5 ≤i≤10 und 20 ≤i≤25, d.h. in diesen Intervallen ist di-

sperse Phase pr¨

asent. Die R¨

ander dieser Intervalle bilden die Phasengrenze Γ. Es soll

nun eine Reihe f′

ikonstruiert werden, deren Phasengrenze Γ′die urspr¨

ungliche Pha-

sengrenze in einem festen Abstand umschließt. Der gew¨

unschte Abstand betrage iδ

Zellbreiten. Zun¨

achst wird durch einmaliges r¨

aumliches Gl¨

atten von fidie Datenreihe

ierzeugt.

i= max[(fi−1+fi+fi+1)/3, fi] (A.19)

0 5 10 15 20 25 30 35

-0,2

0,0

0,2

0,4

0,6

0,8

1,0

f '

' '

Abbildung A.1.: Eindimensionale Verteilung der VoF-Variablen. Kreise: Urspr¨

ungliche

Verteilung fimit der Phasengrenze Γ. ’X’-Symbole: 5-fache Gl¨

attung f5

i. ’+’-Symbole:

Verteilung f′

imit der Phasengrenze Γ′im Abstand von 3 Zellen zu Γ.

120

A.3. Konstruktion einer einh¨

ullenden Fl¨

ache parallel zur Phasengrenze

Durch n-fache Gl¨

attung entsteht

i= max[(fn−1

i−1+fn−1

i+fn−1

i+1 )/3, fn−1

i].(A.20)

Nach n-facher Gl¨

attung enthalten auch Zellen im Abstand nvon der urspr¨

unglichen

Phasengrenze positive Werte. Sei iδ= 3, so wird eine VoF-Verteilung f′

igesucht, deren

Phasengrenzen gegen¨

uber fium jeweils drei Zellen nach außen verschoben sind. Im

Beispiel sei also f′

i>0 in den Intervallen 2 ≤i≤13 und 17 ≤i≤28. Es soll f¨

f′

ialso eine Phasengrenze bei i= 2 liegen. Es muß dann gelten f′

1= 0, f′

3= 1 und

0< f′

2<1. Eine derartiges f′

il¨

aßt sich aus f5

iin einfacher Weise konstruieren:

f′

i=









1 f¨

urf5

i≥fo

0 f¨

urf5

i≤fu

i−fu

fo−fu

f¨

urfu< f5

i< fo,

(A.21)

mit fu:= (f5

1+f5

2)/2 und fo:= (f5

2+f5

3)/2.

Die so konstruierte Verteilung f′

ierf¨

ullt die Anforderung, daß ihre Phasengrenze Γ′

uberall um drei Zellen gegen¨

uber Γ verschoben ist (vgl. Abbildung A.1). Der Algo-

rithmus soll noch einmal f¨

ur einen beliebigen Abstand iδformuliert werden. Sei iγdas

kleinste imit fi>0, so befindet sich eine ”Referenz“-Phasengrenze an der Stelle iγ.

Sei der gesuchte Abstand der neuen Phasengrenze in Einheiten der Zellbreite iδ, so

soll die neue Phasengrenze an der Stelle iγ−iδliegen. Man w¨

ahle ein nmit n > iδund

erh¨

alt durch n-faches r¨

aumliches Gl¨

atten fn

i. Die gesuchte VoF-Verteilung f′

iergibt

sich dann als

f′

i=









1 f¨

urfn

i≥fo

0 f¨

urfn

i≤fu

i−fu

fo−fu

f¨

urfu< fn

i< fo,

(A.22)

mit fu:= (fn

iγ−iδ−1+fn

iγ−iδ)/2 und fo:= (fn

iγ−iδ+fn

iγ−iδ+1)/2.

121

A. Anhang

A.4. Geschwindigkeitsfeld um einen oszillierenden ¨

Oltropfen

In den Abbildungen A.2 bis A.4 ist die zeitliche Entwicklung des Geschwindigkeits-

feldes f¨

ur einen aufsteigenden Mais¨

oltropfen (d=13,2 mm) dargestellt. Zun¨

achst be-

schleunigt der Tropfen geradlinig. W¨

ahrend dieser Phase wird der Nachlauf des Trop-

fens immer l¨

anger, bis er instabil wird. Schließlich l¨

osen sich Wirbel ab, wechselweise

an der linken und der rechten Tropfenseite. Dies geht einher mit der horizontalen

Oszillation. Siehe auch: Kapitel 7.4, Seite 49.

2 4 6 8 10

vof−function

[−]

vmin = 0.207273

vmax = 13.2349

velocity

[cm/s]

t = 1.50020 s

c = 7501

2 4 6 8 10

vof−function

[−]

vmin = 0.0684267

vmax = 13.2401

velocity

[cm/s]

t = 2.50020 s

c = 12501

2 4 6 8 10

vof−function

[−]

vmin = 0.170184

vmax = 13.3201

velocity

[cm/s]

t = 4.00020 s

c = 20001

Abbildung A.2.: Geschwindigkeitsfeld um einen Mais¨

oltropfen (d=13,2 mm) w¨

ahrend

der Beschleunigungsphase. Links: t= 1,5 s. Mitte: t= 2,5 s. Rechts: t= 4,0 s.

122

A.4. Geschwindigkeitsfeld um einen oszillierenden ¨

Oltropfen

2 4 6 8 10

vof−function

[−]

vmin = 0.0574004

vmax = 14.3006

velocity

[cm/s]

t = 6.25000 s

c = 31250

2 4 6 8 10

vof−function

[−]

vmin = 0.0224154

vmax = 14.3597

velocity

[cm/s]

t = 6.45000 s

c = 32250

Abbildung A.3.: Geschwindigkeitsfeld um einen Mais¨

oltropfen (d=13,2 mm) mit peri-

odischer Wirbelabl¨

osung. Links: t= 6,25 s. Rechts: t= 6,45 s.

2 4 6 8 10

vof−function

[−]

vmin = 0.0977254

vmax = 14.1617

velocity

[cm/s]

t = 6.75000 s

c = −31786

2 4 6 8 10

vof−function

[−]

vmin = 0.111000

vmax = 14.1959

velocity

[cm/s]

t = 6.85000 s

c = −31286

Abbildung A.4.: Geschwindigkeitsfeld um einen Mais¨

oltropfen (d=13,2 mm) mit peri-

odischer Wirbelabl¨

osung. Links: t= 6,75 s. Rechts: t= 6,85 s.

123

A. Anhang

124

Literaturverzeichnis

(Alke 2007) Alke, A.: Numerische Modellierung des Einflusses l¨

oslicher Tenside

auf fluide Zweiphasensysteme. Dissertation, Universit¨

at Paderborn 2007.

(Auton 1984) Auton, T. R.: The dynamics of bubbles, drops and particles in mo-

tion in liquids. Ph.D. Thesis, University of Cambridge 1984.

(Auton 1987) Auton, T. R.: The lift force on a spherical body in a rotational flow.

J. Fluid. Mech. 183, S. 199-218 (1987).

(Bothe 2002) Bothe, D.: Computational Fluid Dynamics, Modellierung und An-

wendung. Vorlesungsskript, Universit¨

at Paderborn, Technische Chemie 2002.

(Bothe et al. 2003a) Bothe, D.; Koebe, M.; Warnecke, H.-J.: VoF-Simulation

of the Rise Behavior of Single Air Bubbles with Oxygen Transfer to the Ambient

Liquid. IBW2 Conference on Transport Phenomena with Moving Boundaries. VDI,

S. 1-13 (2003).

(Bothe et al. 2003b) Bothe D.; Koebe, M.; Wielage, K.; Pr¨

uss, J.; Warnecke,

H.-J.: Direct Numerical Simulation of Mass Transfer Between Rising Gas Bubbles

and Water. Sommerfeld; M (Hrsg.): Bubbly Flows - Analysis, Modelling and Cal-

culation. Springer Verlag, Berlin, Heidelberg, New York 2003.

(Bothe et al 2003c) Bothe D.; Koebe, M.; Wielage, K.; Warnecke, H.-J.: VoF-

Simulations of Mass Transfer from Bubbles and Bubble Chains Rising in Aqueous

Solutions. ASME FEDSM 2003, Honolulu, USA, ASME (Veranst.) S. 1-7 (2003)

(Brackbill et al. 1992) Brackbill, J. U.; Kothe, D. B.; Zemach, C.: A Continuum

Method for Modelling Surface Tension. J. Comp. Physics 100, S. 335-354 (1992).

(Brauer 1971) Brauer, H.: Grundlagen der Ein- und Mehrphasenstr¨

omungen. Ver-

lag Sauerl¨

ander, Aarau und Frankfurt/Main 1971.

125

Literaturverzeichnis

(Clift et al. 1978) Clift, R.; Grace, J. R.; Weber, M. E.: Bubbles, Drops and

Particles. Academic Press 1978.

(Drew und Lahey Jr. 1987) Drew, D.; Lahey Jr., R. T.: The virtual mass and

lift force on a sphere in a rotating and straining flow. Int J. Multiphase Flow 13,

S. 113-121 (1987).

(Ervin und Tryggvason 1997) Ervin, E. A.; Tryggvason, G.: The Rise of Bubbles

in Vertical Shear Flow. J. of Fluids Eng. 119, S. 443-449 (1997).

(Fan und Tsuchiya 1990) Fan, L. S.; Tsuchiya, K.: Bubble Wake Dynamics in

Liquids and Liquid-Solid Suspension. Butterworth-Heinemann, Boston 1990.

(Feng und Michaelides 1944) Feng, Z. G.; Michaelides, E. E.: Drag coefficients of

viscous spheres at intermediate and high Reynolds numbers. J. Fluid Eng. - Trans.

ASME, 123(4) S. 841-849 (2001).

(Frohn 1990) Frohn, A.: Numerische Modellierung der Dynamik freier Oberfl¨

achen

in Zweiphasenstr¨

omungen. DFG-Bericht zum Forschungsvorhaben Fr-235/43-1

1999.

(Grossetete 1995) Grossetete C.: Experimental investigation of void profile deve-

lopment in a cylindrical pipe. In Advances in multiphase flow (Serizawa, A.; Fukano,

T.; Bataille, J.; Eds) Amsterdam Elsevier 1995.

(Haberman und Morton 1953) Haberman, W. L.; Morton, R. K.: An experi-

mental investigation of the drag and shape of air bubbles rising in various fluids.

Navy Dept., The David W. Taylor Model Basin 1953.

(Hadamard 1911) Hadamard, J. S.: Mouvement permament lent d’une sph´ere

liquide et visqueuse dans un liquide visqueux. Comptes Rendus des Sc´eances de

l’Acad´emie des Sciences, 152 S. 1735-1738 (1911).

(Haljasmaa et al. 2005) Haljasmaa, I. V.; Vipperman, J. S.; Lynn, R. J.; War-

zinski, R. P.: Control of a fluid particle under simulated deep-ocean conditions in a

high-pressure water tunnel. Rev. Sci. Instrum., 76(2) S. 1-11 (2005).

126

Literaturverzeichnis

(Haljasmaa 2006) Haljasmaa, I.: On the Drag of Fluid and Solid Particles freely

moving in a Continuous Medium. Dissertation, University of Pittsburgh 2006.

(Hirt und Nichols 1981) Hirt, C. W.; Nichols, B. D.: Volume of Fluid (VoF)

method for the dynamics of free Boundaries. J. Comp. Phys. 39, S. 201-225 (1981).

(Ishii 1975) Ishii, M.: Themo-Fluid Dynamic Theory of Two Phase Flow. Eyrolles,

Paris 1975.

(Koebe et al. 2002) Koebe, M.; Bothe, D.; Pr¨

uss, J.; Warnecke, H.-J.: 3D Direct

Numerical Simulation of Air Bubbles in Water at High Reynolds-Numbers. ASME

FEDSM 2002; Montreal; Canada ASME (Veranst.), S. 1-8 (2002).

(Koebe et al. 2003) Koebe, M.; Bothe, D.; Warnecke, H.-J.: Direct Numerical

Simulation of Air Bubbles in Water/Glycerol Mixtures: Shapes and Velocity Fields.

ASME FEDSM 2003, Honolulu, USA, ASME (Veranst.), S. 1-7 (2003).

(Koebe 2004) Koebe, M.: Numerische Simulation aufsteigender Blasen mit und

ohne Stoffaustausch mittels der Volume of Fluid (VoF) Methode. Dissertation, Uni-

versit¨

at Paderborn, 2004.

(Kuipers et al. 2006) Numerical investigation of the lift force acting on single

air bubbles in liquids with different viscosities using a 3D Front Tracking method.

Euromech Colloquium 479, Scheveningen, Niederlande (2006).

(Kurose et al. 2001) Kurose, R.; Misumi, R.; Komori, S.: Drag and lift forces on

a spherical bubble in linear shear flow. Int. J. Multiphase Flow 27, S. 1247-1258

(2001).

(Lafaurie et al. 1994) Lafaurie, B.; Nardonne, C.; Scardovelli, R.; Zaleski, S.; Za-

netti, G.: Modelling, merging and fragmentation in multiphase flow with SURFER.

J. Comp. Phys. 113, S. 134-147 (1994).

(Legrende und Magnaudet 1998) Legrende, D.; Magnaudet, J.: The lift force

on a spherical bubble in a viscous linear shear flow. J. Fluid Mech. 368, S. 81-126

(1998).

127

Literaturverzeichnis

(LeVan und Newman 1976) LeVan, M. D.; Newman, J.: The effect of surfactant

on the terminal and interfacial velocities of a bubble or drop. AIChE J.,22(4), S. 695-

701 (1976).

(Levich 1962) Levich, V. G.: Physicochemical Hydrodynamics. N. J. Enlewood

Cliffs, Prentice Hall 1962.

(Lucas et al. 2005) Lucas, D.; Prasser, H. M.; Manera, A.: Influence of the lift

force on the stability of a bubble column. Chem. Eng. Sci. 60, S. 3609-3619 (2005).

(Maxworthy et al. 1996) Maxworthy, T.; Gnann, C.; K¨

urten, M.; Durst, F.:

Experiments on the Rise of Air Bubbles in Clean Viscous Liquids. Journal of Fluid

Mechanics 321, S. 421-441 (1996).

(Natarajan und Acrivos 1993) Natarajan, R.; Acrivos, A.: The instability of

steady flow past spheres and disks. J. Fluid Mech. 254, S. 323-344 (1993).

(Prasser et al. 2005) Prasser, H. M.; Beyer, M.; Gregor, S.; Lucas, D.; Pietruske,

H.; Sch¨

utz, P.: Weiss, F. P.: Evolution of the Structure of a Gas-Liquid Two-Phase

Flow in a Large Vertical Pipe. NURETH-11, Avinion, France (2005).

(Preukschat 1962) Preukschat, A. W.: Measurements of drag coefficients for fal-

ling and rising spheres in free motion. Calif. Inst. Technol., Pasadena 1962.

(Raymont und Rosant 1999) Raymont, F.; Rosant, J-M.: A numerical and ex-

perimental study of the terminal velocity ans shape of bubbles in viscous liquids.

Chem. Eng. Sci. 55, S. 943-955 (1999).

(Reipschl¨

ager 2002) Reipschl¨

ager, O.: Desintegration von Fluiden im Ultraschall-

feld: Modellierung, Simulation, Experiment. Universit¨

at Paderborn, Dissertation

2002.

(Richard et al. 2007) Richard, A.; Legendre, D.; Magnaudet, J.: Force and torque

on an oblate spheroidal bubble in a linear shear flow. ICMF 2007, Leipzig 2007.

(Rieber und Frohn 1997) Rieber, M.; Frohn, A.: Navier-Stokes simulation of dro-

plet collision dynamics. Proc. 7th Int. Symp. On Comp. Fluid Dynamics, (Zuhang,

F. G. Ed) Bejing, China, S. 520-525 (1997).

128

Literaturverzeichnis

(Rieber 2004) Rieber, M.: Numerische Modellierung der Dynamik freier Grenz-

fl¨

achen in Zweiphasenstr¨

omungen. Fortschritt-Berichte VDI, Reihe 7, Str¨

omungs-

technik 459 Ph.D. thesis, Stuttgart University, Germany 2004.

(Rider und Kothe 1998) Rider, W. J.; Kothe, D. B.: Reconstructing volume

tracking. J. Comp. Phys. 141, S. 112-152 (1998).

(Rusche 2002) Rusche, H.: Computational Fluid Dynamics of Dispersed Two-Phase

Flows at High Phase Fractions. Dissertation, University of London 2002.

(Rybczynski 1911) Rybczynski, W.: On the translatory motion of a fluid sphere in

a viscous medium. Bulletin International de l’Acad´emie Polonaise dews Sciences et

des Lettres, Classe des Sciences Mathe´ematiques et Naturelles, A., S. 40-46 (1911).

(Sabisch 2000) Sabisch, W.: Dreidimensionale numerische Simulation der Dyna-

mik von aufsteigenden Einzelblasen und Blasenschw¨

armen mit einer Volume-of-

Fluid Methode. Forschungszentrum Karlsruhe, Dissertation 2000.

(Sadhal et al. 1997) Sadhal, S. S.; Ayyaswamy, P. S.; Chung, J. N.: Transport

Phenomena with Drops and Bubbles. Springer Verlag, Berlin, Heidelberg, New York

1997.

(Scardovelli und Zaleski 1999) Scardovelli, R.; Zaleski, S.: Direct Numerical

Simulation of Free-Surface Flow and Interfacial Flow. Annual Review of Fluid

Mechanics 31, S. 567-603 (1999).

(Slattery 1999) Slattery, J. C.: Advanced Transport Phenomena. Cambrigde Uni-

versity Press, 1999

(Sousa et al. 2006) Sousa, F. S.; Portela, L. M.; Mudde, R. F.; Mangiavacchi, N.:

Direct Numerical Simulation of Deforable Bubbles in Wall Bounded Schear Flows.

European Conference on Computational Fluid Dynamics, TU Delft, S. 1-13 (2006)

(Stokes 1851) Stokes, G. G.: On the Effect of the Internal Friction of Fluids on

the Motion of Pendulums. Cambridge Philos. Trans., 9 S. 8-106 (1851).

(Tomiyama et al. 1993a) Tomiyama, A.; Sou, A.; Minagawa, H.; Sakaguchi, T.:

Numerical Analysisof a Single Bubble by VoF-Method. JSME International Journal

36, Nr. 1, S. 51-56 (1993).

129

Literaturverzeichnis

(Tomiyama et al. 1995) Tomiyama, A.; Sou, A.; Zun, I; Kanami, A.; Sakaguchi,

T.: Effect of E¨

otv¨

os number and dimensionless liquid volumetric flux on lateral mo-

tion of a bubble in laminar duct flow. In: Advances in Multiphase Flow ( Serizawa,

A.; Fukano, T.; Bataille, J.; Eds) Amsterdam, Elsevier, S. 3-15 (1995).

(Tomiyama et al. 2002) Tomiyama, A.; Tamai, H.; Zun, I.; Hosokawa, S.: Trans-

verse Migration of single Bubbles in Simple Shear Flows. Chem. Eng. Sci. 57(11),

S. 1849-1858 (2002).

(Tomiyama 2004) Tomiyama, A.: Drag, lift and virtual mass forces on a single

bubble. 3th International Symposium on Two-Phase Flow Modeling and Experi-

mentation, Pisa 2004.

(Turton und Levenspiel 1986) Turton, R.; Levenspiel, O.: A short note on the

drag correlation for spheres. Powder Technol. 47, S. 83-86 (1986).

(US Department of Energy Report 1999) US Department of Energy Report:

Carbon sequestation research and developement. DOE/SC/FE-1, available from

NTIS 1999.

(Warzinski et al. 2002) Warzinski, R. P.; Lynn, R. J.; Robertson, A. M.; Haljas-

maa, I. V.: High-pressure water tunnel facility for ocean CO2storage experimenta-

tion. ACS Division of Fuel Chemistry, 47(1), S. 25-26 (2002).

(Wellek et al. 1966) Wellek, R. M.; Agrawal, A. K.; Skelland, A.H.P.: Shapes of

liquid drops moving in liquid media. A.I.Ch.E. Journal, 12, S. 854-862 (1966).

(Zun 1980) Zun, I.: Transverse migration of bubbles influenced by walls in vertical

flow. Int J. Multiphase Flow 6, S. 583-588 (1980).

(Zun 1988) Zun I.: Transition from wall void peaking to core void peaking in turbu-

lent bubbly flow. Proceedings of transient phenomena in multiphase flow. ICHMT

international seminar, Dubrovnik, Croatia, S. 225-233 (1988).

130